diff --git a/bento.info b/bento.info index 99c9453d..9f69677e 100644 --- a/bento.info +++ b/bento.info @@ -93,9 +93,12 @@ Library: Extension: skimage.feature.censure_cy Sources: skimage/feature/censure_cy.pyx - Extension: skimage.feature._brief_cy + Extension: skimage.feature.orb_cy Sources: - skimage/feature/_brief_cy.pyx + skimage/feature/orb_cy.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/doc/examples/plot_brief.py b/doc/examples/plot_brief.py new file mode 100644 index 00000000..47c9ad5f --- /dev/null +++ b/doc/examples/plot_brief.py @@ -0,0 +1,61 @@ +""" +======================= +BRIEF binary descriptor +======================= + +This example demonstrates the BRIEF binary description algorithm. + +The descriptor consists of relatively few bits and can be computed using +a set of intensity difference tests. The short binary descriptor results +in low memory footprint and very efficient matching based on the Hamming +distance metric. + +BRIEF does not provide rotation-invariance. Scale-invariance can be achieved by +detecting and extracting features at different scales. + +""" +from skimage import data +from skimage import transform as tf +from skimage.feature import (match_descriptors, corner_peaks, corner_harris, + plot_matches, BRIEF) +from skimage.color import rgb2gray +import matplotlib.pyplot as plt + + +img1 = rgb2gray(data.lena()) +tform = tf.AffineTransform(scale=(1.2, 1.2), translation=(0, -100)) +img2 = tf.warp(img1, tform) +img3 = tf.rotate(img1, 25) + +keypoints1 = corner_peaks(corner_harris(img1), min_distance=5) +keypoints2 = corner_peaks(corner_harris(img2), min_distance=5) +keypoints3 = corner_peaks(corner_harris(img3), min_distance=5) + +extractor = BRIEF() + +extractor.extract(img1, keypoints1) +keypoints1 = keypoints1[extractor.mask] +descriptors1 = extractor.descriptors + +extractor.extract(img2, keypoints2) +keypoints2 = keypoints2[extractor.mask] +descriptors2 = extractor.descriptors + +extractor.extract(img3, keypoints3) +keypoints3 = keypoints3[extractor.mask] +descriptors3 = extractor.descriptors + +matches12 = match_descriptors(descriptors1, descriptors2, cross_check=True) +matches13 = match_descriptors(descriptors1, descriptors3, cross_check=True) + +fig, ax = plt.subplots(nrows=2, ncols=1) + +plt.gray() + +plot_matches(ax[0], img1, img2, keypoints1, keypoints2, matches12) +ax[0].axis('off') + +plot_matches(ax[1], img1, img3, keypoints1, keypoints3, matches13) +ax[1].axis('off') + +plt.show() diff --git a/doc/examples/plot_censure.py b/doc/examples/plot_censure.py new file mode 100644 index 00000000..c7d70ea5 --- /dev/null +++ b/doc/examples/plot_censure.py @@ -0,0 +1,43 @@ +""" +======================== +CENSURE feature detector +======================== + +The CENSURE feature detector is a scale-invariant center-surround detector +(CENSURE) that claims to outperform other detectors and is capable of real-time +implementation. + +""" +from skimage import data +from skimage import transform as tf +from skimage.feature import CENSURE +from skimage.color import rgb2gray +import matplotlib.pyplot as plt + + +img1 = rgb2gray(data.lena()) +tform = tf.AffineTransform(scale=(1.5, 1.5), rotation=0.5, + translation=(150, -200)) +img2 = tf.warp(img1, tform) + +detector = CENSURE() + +fig, ax = plt.subplots(nrows=1, ncols=2) + +plt.gray() + +detector.detect(img1) + +ax[0].imshow(img1) +ax[0].axis('off') +ax[0].scatter(detector.keypoints[:, 1], detector.keypoints[:, 0], + 2 ** detector.scales, facecolors='none', edgecolors='r') + +detector.detect(img2) + +ax[1].imshow(img2) +ax[1].axis('off') +ax[1].scatter(detector.keypoints[:, 1], detector.keypoints[:, 0], + 2 ** detector.scales, facecolors='none', edgecolors='r') + +plt.show() diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py index bae6d6e2..5812fa37 100644 --- a/doc/examples/plot_matching.py +++ b/doc/examples/plot_matching.py @@ -27,7 +27,8 @@ from matplotlib import pyplot as plt from skimage import data from skimage.util import img_as_float -from skimage.feature import corner_harris, corner_subpix, corner_peaks +from skimage.feature import (corner_harris, corner_subpix, corner_peaks, + plot_matches) from skimage.transform import warp, AffineTransform from skimage.exposure import rescale_intensity from skimage.color import rgb2gray @@ -117,28 +118,21 @@ print(tform.scale, tform.translation, tform.rotation) print(model.scale, model.translation, model.rotation) print(model_robust.scale, model_robust.translation, model_robust.rotation) - -# visualize correspondences -img_combined = np.concatenate((img_orig_gray, img_warped_gray), axis=1) - +# visualize correspondence fig, ax = plt.subplots(nrows=2, ncols=1) + plt.gray() -ax[0].imshow(img_combined, interpolation='nearest') +inlier_idxs = np.nonzero(inliers)[0] +plot_matches(ax[0], img_orig_gray, img_warped_gray, src, dst, + np.column_stack((inlier_idxs, inlier_idxs)), matches_color='b') ax[0].axis('off') -ax[0].axis((0, 400, 200, 0)) ax[0].set_title('Correct correspondences') -ax[1].imshow(img_combined, interpolation='nearest') + +outlier_idxs = np.nonzero(outliers)[0] +plot_matches(ax[1], img_orig_gray, img_warped_gray, src, dst, + np.column_stack((outlier_idxs, outlier_idxs)), matches_color='r') ax[1].axis('off') -ax[1].axis((0, 400, 200, 0)) ax[1].set_title('Faulty correspondences') - -for ax_idx, (m, color) in enumerate(((inliers, 'g'), (outliers, 'r'))): - ax[ax_idx].plot((src[m, 1], dst[m, 1] + 200), (src[m, 0], dst[m, 0]), '-', - color=color) - ax[ax_idx].plot(src[m, 1], src[m, 0], '.', markersize=10, color=color) - ax[ax_idx].plot(dst[m, 1] + 200, dst[m, 0], '.', markersize=10, - color=color) - plt.show() diff --git a/doc/examples/plot_orb.py b/doc/examples/plot_orb.py new file mode 100644 index 00000000..1a73fc7f --- /dev/null +++ b/doc/examples/plot_orb.py @@ -0,0 +1,56 @@ +""" +========================================== +ORB feature detector and binary descriptor +========================================== + +This example demonstrates the ORB feature detection and binary description +algorithm. It uses an oriented FAST detection method and the rotated BRIEF +descriptors. + +Unlike BRIEF, ORB is comparatively scale- and rotation-invariant while still +employing the very efficient Hamming distance metric for matching. As such, it +is preferred for real-time applications. + +""" +from skimage import data +from skimage import transform as tf +from skimage.feature import (match_descriptors, corner_harris, + corner_peaks, ORB, plot_matches) +from skimage.color import rgb2gray +import matplotlib.pyplot as plt + + +img1 = rgb2gray(data.lena()) +img2 = tf.rotate(img1, 180) +tform = tf.AffineTransform(scale=(1.3, 1.1), rotation=0.5, + translation=(0, -200)) +img3 = tf.warp(img1, tform) + +descriptor_extractor = ORB(n_keypoints=200) + +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, cross_check=True) +matches13 = match_descriptors(descriptors1, descriptors3, cross_check=True) + +fig, ax = plt.subplots(nrows=2, ncols=1) + +plt.gray() + +plot_matches(ax[0], img1, img2, keypoints1, keypoints2, matches12) +ax[0].axis('off') + +plot_matches(ax[1], img1, img3, keypoints1, keypoints3, matches13) +ax[1].axis('off') + +plt.show() diff --git a/skimage/data/orb_descriptor_positions.txt b/skimage/data/orb_descriptor_positions.txt new file mode 100644 index 00000000..541f972f --- /dev/null +++ b/skimage/data/orb_descriptor_positions.txt @@ -0,0 +1,256 @@ +8.000000000000000000e+00 -3.000000000000000000e+00 9.000000000000000000e+00 5.000000000000000000e+00 +4.000000000000000000e+00 2.000000000000000000e+00 7.000000000000000000e+00 -1.200000000000000000e+01 +-1.100000000000000000e+01 9.000000000000000000e+00 -8.000000000000000000e+00 2.000000000000000000e+00 +7.000000000000000000e+00 -1.200000000000000000e+01 1.200000000000000000e+01 -1.300000000000000000e+01 +2.000000000000000000e+00 -1.300000000000000000e+01 2.000000000000000000e+00 1.200000000000000000e+01 +1.000000000000000000e+00 -7.000000000000000000e+00 1.000000000000000000e+00 6.000000000000000000e+00 +-2.000000000000000000e+00 -1.000000000000000000e+01 -2.000000000000000000e+00 -4.000000000000000000e+00 +-1.300000000000000000e+01 -1.300000000000000000e+01 -1.100000000000000000e+01 -8.000000000000000000e+00 +-1.300000000000000000e+01 -3.000000000000000000e+00 -1.200000000000000000e+01 -9.000000000000000000e+00 +1.000000000000000000e+01 4.000000000000000000e+00 1.100000000000000000e+01 9.000000000000000000e+00 +-1.300000000000000000e+01 -8.000000000000000000e+00 -8.000000000000000000e+00 -9.000000000000000000e+00 +-1.100000000000000000e+01 7.000000000000000000e+00 -9.000000000000000000e+00 1.200000000000000000e+01 +7.000000000000000000e+00 7.000000000000000000e+00 1.200000000000000000e+01 6.000000000000000000e+00 +-4.000000000000000000e+00 -5.000000000000000000e+00 -3.000000000000000000e+00 0.000000000000000000e+00 +-1.300000000000000000e+01 2.000000000000000000e+00 -1.200000000000000000e+01 -3.000000000000000000e+00 +-9.000000000000000000e+00 0.000000000000000000e+00 -7.000000000000000000e+00 5.000000000000000000e+00 +1.200000000000000000e+01 -6.000000000000000000e+00 1.200000000000000000e+01 -1.000000000000000000e+00 +-3.000000000000000000e+00 6.000000000000000000e+00 -2.000000000000000000e+00 1.200000000000000000e+01 +-6.000000000000000000e+00 -1.300000000000000000e+01 -4.000000000000000000e+00 -8.000000000000000000e+00 +1.100000000000000000e+01 -1.300000000000000000e+01 1.200000000000000000e+01 -8.000000000000000000e+00 +4.000000000000000000e+00 7.000000000000000000e+00 5.000000000000000000e+00 1.000000000000000000e+00 +5.000000000000000000e+00 -3.000000000000000000e+00 1.000000000000000000e+01 -3.000000000000000000e+00 +3.000000000000000000e+00 -7.000000000000000000e+00 6.000000000000000000e+00 1.200000000000000000e+01 +-8.000000000000000000e+00 -7.000000000000000000e+00 -6.000000000000000000e+00 -2.000000000000000000e+00 +-2.000000000000000000e+00 1.100000000000000000e+01 -1.000000000000000000e+00 -1.000000000000000000e+01 +-1.300000000000000000e+01 1.200000000000000000e+01 -8.000000000000000000e+00 1.000000000000000000e+01 +-7.000000000000000000e+00 3.000000000000000000e+00 -5.000000000000000000e+00 -3.000000000000000000e+00 +-4.000000000000000000e+00 2.000000000000000000e+00 -3.000000000000000000e+00 7.000000000000000000e+00 +-1.000000000000000000e+01 -1.200000000000000000e+01 -6.000000000000000000e+00 1.100000000000000000e+01 +5.000000000000000000e+00 -1.200000000000000000e+01 6.000000000000000000e+00 -7.000000000000000000e+00 +5.000000000000000000e+00 -6.000000000000000000e+00 7.000000000000000000e+00 -1.000000000000000000e+00 +1.000000000000000000e+00 0.000000000000000000e+00 4.000000000000000000e+00 -5.000000000000000000e+00 +9.000000000000000000e+00 1.100000000000000000e+01 1.100000000000000000e+01 -1.300000000000000000e+01 +4.000000000000000000e+00 7.000000000000000000e+00 4.000000000000000000e+00 1.200000000000000000e+01 +2.000000000000000000e+00 -1.000000000000000000e+00 4.000000000000000000e+00 4.000000000000000000e+00 +-4.000000000000000000e+00 -1.200000000000000000e+01 -2.000000000000000000e+00 7.000000000000000000e+00 +-8.000000000000000000e+00 -5.000000000000000000e+00 -7.000000000000000000e+00 -1.000000000000000000e+01 +4.000000000000000000e+00 1.100000000000000000e+01 9.000000000000000000e+00 1.200000000000000000e+01 +0.000000000000000000e+00 -8.000000000000000000e+00 1.000000000000000000e+00 -1.300000000000000000e+01 +-1.300000000000000000e+01 -2.000000000000000000e+00 -8.000000000000000000e+00 2.000000000000000000e+00 +-3.000000000000000000e+00 -2.000000000000000000e+00 -2.000000000000000000e+00 3.000000000000000000e+00 +-6.000000000000000000e+00 9.000000000000000000e+00 -4.000000000000000000e+00 -9.000000000000000000e+00 +8.000000000000000000e+00 1.200000000000000000e+01 1.000000000000000000e+01 7.000000000000000000e+00 +0.000000000000000000e+00 9.000000000000000000e+00 1.000000000000000000e+00 3.000000000000000000e+00 +7.000000000000000000e+00 -5.000000000000000000e+00 1.100000000000000000e+01 -1.000000000000000000e+01 +-1.300000000000000000e+01 -6.000000000000000000e+00 -1.100000000000000000e+01 0.000000000000000000e+00 +1.000000000000000000e+01 7.000000000000000000e+00 1.200000000000000000e+01 1.000000000000000000e+00 +-6.000000000000000000e+00 -3.000000000000000000e+00 -6.000000000000000000e+00 1.200000000000000000e+01 +1.000000000000000000e+01 -9.000000000000000000e+00 1.200000000000000000e+01 -4.000000000000000000e+00 +-1.300000000000000000e+01 8.000000000000000000e+00 -8.000000000000000000e+00 -1.200000000000000000e+01 +-1.300000000000000000e+01 0.000000000000000000e+00 -8.000000000000000000e+00 -4.000000000000000000e+00 +3.000000000000000000e+00 3.000000000000000000e+00 7.000000000000000000e+00 8.000000000000000000e+00 +5.000000000000000000e+00 7.000000000000000000e+00 1.000000000000000000e+01 -7.000000000000000000e+00 +-1.000000000000000000e+00 7.000000000000000000e+00 1.000000000000000000e+00 -1.200000000000000000e+01 +3.000000000000000000e+00 -1.000000000000000000e+01 5.000000000000000000e+00 6.000000000000000000e+00 +2.000000000000000000e+00 -4.000000000000000000e+00 3.000000000000000000e+00 -1.000000000000000000e+01 +-1.300000000000000000e+01 0.000000000000000000e+00 -1.300000000000000000e+01 5.000000000000000000e+00 +-1.300000000000000000e+01 -7.000000000000000000e+00 -1.200000000000000000e+01 1.200000000000000000e+01 +-1.300000000000000000e+01 3.000000000000000000e+00 -1.100000000000000000e+01 8.000000000000000000e+00 +-7.000000000000000000e+00 1.200000000000000000e+01 -4.000000000000000000e+00 7.000000000000000000e+00 +6.000000000000000000e+00 -1.000000000000000000e+01 1.200000000000000000e+01 8.000000000000000000e+00 +-9.000000000000000000e+00 -1.000000000000000000e+00 -7.000000000000000000e+00 -6.000000000000000000e+00 +-2.000000000000000000e+00 -5.000000000000000000e+00 0.000000000000000000e+00 1.200000000000000000e+01 +-1.200000000000000000e+01 5.000000000000000000e+00 -7.000000000000000000e+00 5.000000000000000000e+00 +3.000000000000000000e+00 -1.000000000000000000e+01 8.000000000000000000e+00 -1.300000000000000000e+01 +-7.000000000000000000e+00 -7.000000000000000000e+00 -4.000000000000000000e+00 5.000000000000000000e+00 +-3.000000000000000000e+00 -2.000000000000000000e+00 -1.000000000000000000e+00 -7.000000000000000000e+00 +2.000000000000000000e+00 9.000000000000000000e+00 5.000000000000000000e+00 -1.100000000000000000e+01 +-1.100000000000000000e+01 -1.300000000000000000e+01 -5.000000000000000000e+00 -1.300000000000000000e+01 +-1.000000000000000000e+00 6.000000000000000000e+00 0.000000000000000000e+00 -1.000000000000000000e+00 +5.000000000000000000e+00 -3.000000000000000000e+00 5.000000000000000000e+00 2.000000000000000000e+00 +-4.000000000000000000e+00 -1.300000000000000000e+01 -4.000000000000000000e+00 1.200000000000000000e+01 +-9.000000000000000000e+00 -6.000000000000000000e+00 -9.000000000000000000e+00 6.000000000000000000e+00 +-1.200000000000000000e+01 -1.000000000000000000e+01 -8.000000000000000000e+00 -4.000000000000000000e+00 +1.000000000000000000e+01 2.000000000000000000e+00 1.200000000000000000e+01 -3.000000000000000000e+00 +7.000000000000000000e+00 1.200000000000000000e+01 1.200000000000000000e+01 1.200000000000000000e+01 +-7.000000000000000000e+00 -1.300000000000000000e+01 -6.000000000000000000e+00 5.000000000000000000e+00 +-4.000000000000000000e+00 9.000000000000000000e+00 -3.000000000000000000e+00 4.000000000000000000e+00 +7.000000000000000000e+00 -1.000000000000000000e+00 1.200000000000000000e+01 2.000000000000000000e+00 +-7.000000000000000000e+00 6.000000000000000000e+00 -5.000000000000000000e+00 1.000000000000000000e+00 +-1.300000000000000000e+01 1.100000000000000000e+01 -1.200000000000000000e+01 5.000000000000000000e+00 +-3.000000000000000000e+00 7.000000000000000000e+00 -2.000000000000000000e+00 -6.000000000000000000e+00 +7.000000000000000000e+00 -8.000000000000000000e+00 1.200000000000000000e+01 -7.000000000000000000e+00 +-1.300000000000000000e+01 -7.000000000000000000e+00 -1.100000000000000000e+01 -1.200000000000000000e+01 +1.000000000000000000e+00 -3.000000000000000000e+00 1.200000000000000000e+01 1.200000000000000000e+01 +2.000000000000000000e+00 -6.000000000000000000e+00 3.000000000000000000e+00 0.000000000000000000e+00 +-4.000000000000000000e+00 3.000000000000000000e+00 -2.000000000000000000e+00 -1.300000000000000000e+01 +-1.000000000000000000e+00 -1.300000000000000000e+01 1.000000000000000000e+00 9.000000000000000000e+00 +7.000000000000000000e+00 1.000000000000000000e+00 8.000000000000000000e+00 -6.000000000000000000e+00 +1.000000000000000000e+00 -1.000000000000000000e+00 3.000000000000000000e+00 1.200000000000000000e+01 +9.000000000000000000e+00 1.000000000000000000e+00 1.200000000000000000e+01 6.000000000000000000e+00 +-1.000000000000000000e+00 -9.000000000000000000e+00 -1.000000000000000000e+00 3.000000000000000000e+00 +-1.300000000000000000e+01 -1.300000000000000000e+01 -1.000000000000000000e+01 5.000000000000000000e+00 +7.000000000000000000e+00 7.000000000000000000e+00 1.000000000000000000e+01 1.200000000000000000e+01 +1.200000000000000000e+01 -5.000000000000000000e+00 1.200000000000000000e+01 9.000000000000000000e+00 +6.000000000000000000e+00 3.000000000000000000e+00 7.000000000000000000e+00 1.100000000000000000e+01 +5.000000000000000000e+00 -1.300000000000000000e+01 6.000000000000000000e+00 1.000000000000000000e+01 +2.000000000000000000e+00 -1.200000000000000000e+01 2.000000000000000000e+00 3.000000000000000000e+00 +3.000000000000000000e+00 8.000000000000000000e+00 4.000000000000000000e+00 -6.000000000000000000e+00 +2.000000000000000000e+00 6.000000000000000000e+00 1.200000000000000000e+01 -1.300000000000000000e+01 +9.000000000000000000e+00 -1.200000000000000000e+01 1.000000000000000000e+01 3.000000000000000000e+00 +-8.000000000000000000e+00 4.000000000000000000e+00 -7.000000000000000000e+00 9.000000000000000000e+00 +-1.100000000000000000e+01 1.200000000000000000e+01 -4.000000000000000000e+00 -6.000000000000000000e+00 +1.000000000000000000e+00 1.200000000000000000e+01 2.000000000000000000e+00 -8.000000000000000000e+00 +6.000000000000000000e+00 -9.000000000000000000e+00 7.000000000000000000e+00 -4.000000000000000000e+00 +2.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 -2.000000000000000000e+00 +6.000000000000000000e+00 3.000000000000000000e+00 1.100000000000000000e+01 0.000000000000000000e+00 +3.000000000000000000e+00 -3.000000000000000000e+00 8.000000000000000000e+00 -8.000000000000000000e+00 +7.000000000000000000e+00 8.000000000000000000e+00 9.000000000000000000e+00 3.000000000000000000e+00 +-1.100000000000000000e+01 -5.000000000000000000e+00 -6.000000000000000000e+00 -4.000000000000000000e+00 +-1.000000000000000000e+01 1.100000000000000000e+01 -5.000000000000000000e+00 1.000000000000000000e+01 +-5.000000000000000000e+00 -8.000000000000000000e+00 -3.000000000000000000e+00 1.200000000000000000e+01 +-1.000000000000000000e+01 5.000000000000000000e+00 -9.000000000000000000e+00 0.000000000000000000e+00 +8.000000000000000000e+00 -1.000000000000000000e+00 1.200000000000000000e+01 -6.000000000000000000e+00 +4.000000000000000000e+00 -6.000000000000000000e+00 6.000000000000000000e+00 -1.100000000000000000e+01 +-1.000000000000000000e+01 1.200000000000000000e+01 -8.000000000000000000e+00 7.000000000000000000e+00 +4.000000000000000000e+00 -2.000000000000000000e+00 6.000000000000000000e+00 7.000000000000000000e+00 +-2.000000000000000000e+00 0.000000000000000000e+00 -2.000000000000000000e+00 1.200000000000000000e+01 +-5.000000000000000000e+00 -8.000000000000000000e+00 -5.000000000000000000e+00 2.000000000000000000e+00 +7.000000000000000000e+00 -6.000000000000000000e+00 1.000000000000000000e+01 1.200000000000000000e+01 +-9.000000000000000000e+00 -1.300000000000000000e+01 -8.000000000000000000e+00 -8.000000000000000000e+00 +-5.000000000000000000e+00 -1.300000000000000000e+01 -5.000000000000000000e+00 -2.000000000000000000e+00 +8.000000000000000000e+00 -8.000000000000000000e+00 9.000000000000000000e+00 -1.300000000000000000e+01 +-9.000000000000000000e+00 -1.100000000000000000e+01 -9.000000000000000000e+00 0.000000000000000000e+00 +1.000000000000000000e+00 -8.000000000000000000e+00 1.000000000000000000e+00 -2.000000000000000000e+00 +7.000000000000000000e+00 -4.000000000000000000e+00 9.000000000000000000e+00 1.000000000000000000e+00 +-2.000000000000000000e+00 1.000000000000000000e+00 -1.000000000000000000e+00 -4.000000000000000000e+00 +1.100000000000000000e+01 -6.000000000000000000e+00 1.200000000000000000e+01 -1.100000000000000000e+01 +-1.200000000000000000e+01 -9.000000000000000000e+00 -6.000000000000000000e+00 4.000000000000000000e+00 +3.000000000000000000e+00 7.000000000000000000e+00 7.000000000000000000e+00 1.200000000000000000e+01 +5.000000000000000000e+00 5.000000000000000000e+00 1.000000000000000000e+01 8.000000000000000000e+00 +0.000000000000000000e+00 -4.000000000000000000e+00 2.000000000000000000e+00 8.000000000000000000e+00 +-9.000000000000000000e+00 1.200000000000000000e+01 -5.000000000000000000e+00 -1.300000000000000000e+01 +0.000000000000000000e+00 7.000000000000000000e+00 2.000000000000000000e+00 1.200000000000000000e+01 +-1.000000000000000000e+00 2.000000000000000000e+00 1.000000000000000000e+00 7.000000000000000000e+00 +5.000000000000000000e+00 1.100000000000000000e+01 7.000000000000000000e+00 -9.000000000000000000e+00 +3.000000000000000000e+00 5.000000000000000000e+00 6.000000000000000000e+00 -8.000000000000000000e+00 +-1.300000000000000000e+01 -4.000000000000000000e+00 -8.000000000000000000e+00 9.000000000000000000e+00 +-5.000000000000000000e+00 9.000000000000000000e+00 -3.000000000000000000e+00 -3.000000000000000000e+00 +-4.000000000000000000e+00 -7.000000000000000000e+00 -3.000000000000000000e+00 -1.200000000000000000e+01 +6.000000000000000000e+00 5.000000000000000000e+00 8.000000000000000000e+00 0.000000000000000000e+00 +-7.000000000000000000e+00 6.000000000000000000e+00 -6.000000000000000000e+00 1.200000000000000000e+01 +-1.300000000000000000e+01 6.000000000000000000e+00 -5.000000000000000000e+00 -2.000000000000000000e+00 +1.000000000000000000e+00 -1.000000000000000000e+01 3.000000000000000000e+00 1.000000000000000000e+01 +4.000000000000000000e+00 1.000000000000000000e+00 8.000000000000000000e+00 -4.000000000000000000e+00 +-2.000000000000000000e+00 -2.000000000000000000e+00 2.000000000000000000e+00 -1.300000000000000000e+01 +2.000000000000000000e+00 -1.200000000000000000e+01 1.200000000000000000e+01 1.200000000000000000e+01 +-2.000000000000000000e+00 -1.300000000000000000e+01 0.000000000000000000e+00 -6.000000000000000000e+00 +4.000000000000000000e+00 1.000000000000000000e+00 9.000000000000000000e+00 3.000000000000000000e+00 +-6.000000000000000000e+00 -1.000000000000000000e+01 -3.000000000000000000e+00 -5.000000000000000000e+00 +-3.000000000000000000e+00 -1.300000000000000000e+01 -1.000000000000000000e+00 1.000000000000000000e+00 +7.000000000000000000e+00 5.000000000000000000e+00 1.200000000000000000e+01 -1.100000000000000000e+01 +4.000000000000000000e+00 -2.000000000000000000e+00 5.000000000000000000e+00 -7.000000000000000000e+00 +-1.300000000000000000e+01 9.000000000000000000e+00 -9.000000000000000000e+00 -5.000000000000000000e+00 +7.000000000000000000e+00 1.000000000000000000e+00 8.000000000000000000e+00 6.000000000000000000e+00 +7.000000000000000000e+00 -8.000000000000000000e+00 7.000000000000000000e+00 6.000000000000000000e+00 +-7.000000000000000000e+00 -4.000000000000000000e+00 -7.000000000000000000e+00 1.000000000000000000e+00 +-8.000000000000000000e+00 1.100000000000000000e+01 -7.000000000000000000e+00 -8.000000000000000000e+00 +-1.300000000000000000e+01 6.000000000000000000e+00 -1.200000000000000000e+01 -8.000000000000000000e+00 +2.000000000000000000e+00 4.000000000000000000e+00 3.000000000000000000e+00 9.000000000000000000e+00 +1.000000000000000000e+01 -5.000000000000000000e+00 1.200000000000000000e+01 3.000000000000000000e+00 +-6.000000000000000000e+00 -5.000000000000000000e+00 -6.000000000000000000e+00 7.000000000000000000e+00 +8.000000000000000000e+00 -3.000000000000000000e+00 9.000000000000000000e+00 -8.000000000000000000e+00 +2.000000000000000000e+00 -1.200000000000000000e+01 2.000000000000000000e+00 8.000000000000000000e+00 +-1.100000000000000000e+01 -2.000000000000000000e+00 -1.000000000000000000e+01 3.000000000000000000e+00 +-1.200000000000000000e+01 -1.300000000000000000e+01 -7.000000000000000000e+00 -9.000000000000000000e+00 +-1.100000000000000000e+01 0.000000000000000000e+00 -1.000000000000000000e+01 -5.000000000000000000e+00 +5.000000000000000000e+00 -3.000000000000000000e+00 1.100000000000000000e+01 8.000000000000000000e+00 +-2.000000000000000000e+00 -1.300000000000000000e+01 -1.000000000000000000e+00 1.200000000000000000e+01 +-1.000000000000000000e+00 -8.000000000000000000e+00 0.000000000000000000e+00 9.000000000000000000e+00 +-1.300000000000000000e+01 -1.100000000000000000e+01 -1.200000000000000000e+01 -5.000000000000000000e+00 +-1.000000000000000000e+01 -2.000000000000000000e+00 -1.000000000000000000e+01 1.100000000000000000e+01 +-3.000000000000000000e+00 9.000000000000000000e+00 -2.000000000000000000e+00 -1.300000000000000000e+01 +2.000000000000000000e+00 -3.000000000000000000e+00 3.000000000000000000e+00 2.000000000000000000e+00 +-9.000000000000000000e+00 -1.300000000000000000e+01 -4.000000000000000000e+00 0.000000000000000000e+00 +-4.000000000000000000e+00 6.000000000000000000e+00 -3.000000000000000000e+00 -1.000000000000000000e+01 +-4.000000000000000000e+00 1.200000000000000000e+01 -2.000000000000000000e+00 -7.000000000000000000e+00 +-6.000000000000000000e+00 -1.100000000000000000e+01 -4.000000000000000000e+00 9.000000000000000000e+00 +6.000000000000000000e+00 -3.000000000000000000e+00 6.000000000000000000e+00 1.100000000000000000e+01 +-1.300000000000000000e+01 1.100000000000000000e+01 -5.000000000000000000e+00 5.000000000000000000e+00 +1.100000000000000000e+01 1.100000000000000000e+01 1.200000000000000000e+01 6.000000000000000000e+00 +7.000000000000000000e+00 -5.000000000000000000e+00 1.200000000000000000e+01 -2.000000000000000000e+00 +-1.000000000000000000e+00 1.200000000000000000e+01 0.000000000000000000e+00 7.000000000000000000e+00 +-4.000000000000000000e+00 -8.000000000000000000e+00 -3.000000000000000000e+00 -2.000000000000000000e+00 +-7.000000000000000000e+00 1.000000000000000000e+00 -6.000000000000000000e+00 7.000000000000000000e+00 +-1.300000000000000000e+01 -1.200000000000000000e+01 -8.000000000000000000e+00 -1.300000000000000000e+01 +-7.000000000000000000e+00 -2.000000000000000000e+00 -6.000000000000000000e+00 -8.000000000000000000e+00 +-8.000000000000000000e+00 5.000000000000000000e+00 -6.000000000000000000e+00 -9.000000000000000000e+00 +-5.000000000000000000e+00 -1.000000000000000000e+00 -4.000000000000000000e+00 5.000000000000000000e+00 +-1.300000000000000000e+01 7.000000000000000000e+00 -8.000000000000000000e+00 1.000000000000000000e+01 +1.000000000000000000e+00 5.000000000000000000e+00 5.000000000000000000e+00 -1.300000000000000000e+01 +1.000000000000000000e+00 0.000000000000000000e+00 1.000000000000000000e+01 -1.300000000000000000e+01 +9.000000000000000000e+00 1.200000000000000000e+01 1.000000000000000000e+01 -1.000000000000000000e+00 +5.000000000000000000e+00 -8.000000000000000000e+00 1.000000000000000000e+01 -9.000000000000000000e+00 +-1.000000000000000000e+00 1.100000000000000000e+01 1.000000000000000000e+00 -1.300000000000000000e+01 +-9.000000000000000000e+00 -3.000000000000000000e+00 -6.000000000000000000e+00 2.000000000000000000e+00 +-1.000000000000000000e+00 -1.000000000000000000e+01 1.000000000000000000e+00 1.200000000000000000e+01 +-1.300000000000000000e+01 1.000000000000000000e+00 -8.000000000000000000e+00 -1.000000000000000000e+01 +8.000000000000000000e+00 -1.100000000000000000e+01 1.000000000000000000e+01 -6.000000000000000000e+00 +2.000000000000000000e+00 -1.300000000000000000e+01 3.000000000000000000e+00 -6.000000000000000000e+00 +7.000000000000000000e+00 -1.300000000000000000e+01 1.200000000000000000e+01 -9.000000000000000000e+00 +-1.000000000000000000e+01 -1.000000000000000000e+01 -5.000000000000000000e+00 -7.000000000000000000e+00 +-1.000000000000000000e+01 -8.000000000000000000e+00 -8.000000000000000000e+00 -1.300000000000000000e+01 +4.000000000000000000e+00 -6.000000000000000000e+00 8.000000000000000000e+00 5.000000000000000000e+00 +3.000000000000000000e+00 1.200000000000000000e+01 8.000000000000000000e+00 -1.300000000000000000e+01 +-4.000000000000000000e+00 2.000000000000000000e+00 -3.000000000000000000e+00 -3.000000000000000000e+00 +5.000000000000000000e+00 -1.300000000000000000e+01 1.000000000000000000e+01 -1.200000000000000000e+01 +4.000000000000000000e+00 -1.300000000000000000e+01 5.000000000000000000e+00 -1.000000000000000000e+00 +-9.000000000000000000e+00 9.000000000000000000e+00 -4.000000000000000000e+00 3.000000000000000000e+00 +0.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 -9.000000000000000000e+00 +-1.200000000000000000e+01 1.000000000000000000e+00 -6.000000000000000000e+00 1.000000000000000000e+00 +3.000000000000000000e+00 2.000000000000000000e+00 4.000000000000000000e+00 -8.000000000000000000e+00 +-1.000000000000000000e+01 -1.000000000000000000e+01 -1.000000000000000000e+01 9.000000000000000000e+00 +8.000000000000000000e+00 -1.300000000000000000e+01 1.200000000000000000e+01 1.200000000000000000e+01 +-8.000000000000000000e+00 -1.200000000000000000e+01 -6.000000000000000000e+00 -5.000000000000000000e+00 +2.000000000000000000e+00 2.000000000000000000e+00 3.000000000000000000e+00 7.000000000000000000e+00 +1.000000000000000000e+01 6.000000000000000000e+00 1.100000000000000000e+01 -8.000000000000000000e+00 +6.000000000000000000e+00 8.000000000000000000e+00 8.000000000000000000e+00 -1.200000000000000000e+01 +-7.000000000000000000e+00 1.000000000000000000e+01 -6.000000000000000000e+00 5.000000000000000000e+00 +-3.000000000000000000e+00 -9.000000000000000000e+00 -3.000000000000000000e+00 9.000000000000000000e+00 +-1.000000000000000000e+00 -1.300000000000000000e+01 -1.000000000000000000e+00 5.000000000000000000e+00 +-3.000000000000000000e+00 -7.000000000000000000e+00 -3.000000000000000000e+00 4.000000000000000000e+00 +-8.000000000000000000e+00 -2.000000000000000000e+00 -8.000000000000000000e+00 3.000000000000000000e+00 +4.000000000000000000e+00 2.000000000000000000e+00 1.200000000000000000e+01 1.200000000000000000e+01 +2.000000000000000000e+00 -5.000000000000000000e+00 3.000000000000000000e+00 1.100000000000000000e+01 +6.000000000000000000e+00 -9.000000000000000000e+00 1.100000000000000000e+01 -1.300000000000000000e+01 +3.000000000000000000e+00 -1.000000000000000000e+00 7.000000000000000000e+00 1.200000000000000000e+01 +1.100000000000000000e+01 -1.000000000000000000e+00 1.200000000000000000e+01 4.000000000000000000e+00 +-3.000000000000000000e+00 0.000000000000000000e+00 -3.000000000000000000e+00 6.000000000000000000e+00 +4.000000000000000000e+00 -1.100000000000000000e+01 4.000000000000000000e+00 1.200000000000000000e+01 +2.000000000000000000e+00 -4.000000000000000000e+00 2.000000000000000000e+00 1.000000000000000000e+00 +-1.000000000000000000e+01 -6.000000000000000000e+00 -8.000000000000000000e+00 1.000000000000000000e+00 +-1.300000000000000000e+01 7.000000000000000000e+00 -1.100000000000000000e+01 1.000000000000000000e+00 +-1.300000000000000000e+01 1.200000000000000000e+01 -1.100000000000000000e+01 -1.300000000000000000e+01 +6.000000000000000000e+00 0.000000000000000000e+00 1.100000000000000000e+01 -1.300000000000000000e+01 +0.000000000000000000e+00 -1.000000000000000000e+00 1.000000000000000000e+00 4.000000000000000000e+00 +-1.300000000000000000e+01 3.000000000000000000e+00 -9.000000000000000000e+00 -2.000000000000000000e+00 +-9.000000000000000000e+00 8.000000000000000000e+00 -6.000000000000000000e+00 -3.000000000000000000e+00 +-1.300000000000000000e+01 -6.000000000000000000e+00 -8.000000000000000000e+00 -2.000000000000000000e+00 +5.000000000000000000e+00 -9.000000000000000000e+00 8.000000000000000000e+00 1.000000000000000000e+01 +2.000000000000000000e+00 7.000000000000000000e+00 3.000000000000000000e+00 -9.000000000000000000e+00 +-1.000000000000000000e+00 -6.000000000000000000e+00 -1.000000000000000000e+00 -1.000000000000000000e+00 +9.000000000000000000e+00 5.000000000000000000e+00 1.100000000000000000e+01 -2.000000000000000000e+00 +1.100000000000000000e+01 -3.000000000000000000e+00 1.200000000000000000e+01 -8.000000000000000000e+00 +3.000000000000000000e+00 0.000000000000000000e+00 3.000000000000000000e+00 5.000000000000000000e+00 +-1.000000000000000000e+00 4.000000000000000000e+00 0.000000000000000000e+00 1.000000000000000000e+01 +3.000000000000000000e+00 -6.000000000000000000e+00 4.000000000000000000e+00 5.000000000000000000e+00 +-1.300000000000000000e+01 0.000000000000000000e+00 -1.000000000000000000e+01 5.000000000000000000e+00 +5.000000000000000000e+00 8.000000000000000000e+00 1.200000000000000000e+01 1.100000000000000000e+01 +8.000000000000000000e+00 9.000000000000000000e+00 9.000000000000000000e+00 -6.000000000000000000e+00 +7.000000000000000000e+00 -4.000000000000000000e+00 8.000000000000000000e+00 -1.200000000000000000e+01 +-1.000000000000000000e+01 4.000000000000000000e+00 -1.000000000000000000e+01 9.000000000000000000e+00 +7.000000000000000000e+00 3.000000000000000000e+00 1.200000000000000000e+01 4.000000000000000000e+00 +9.000000000000000000e+00 -7.000000000000000000e+00 1.000000000000000000e+01 -2.000000000000000000e+00 +7.000000000000000000e+00 0.000000000000000000e+00 1.200000000000000000e+01 -2.000000000000000000e+00 +-1.000000000000000000e+00 -6.000000000000000000e+00 0.000000000000000000e+00 -1.100000000000000000e+01 diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 4a6518d6..1b1cb189 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -4,9 +4,16 @@ 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_cy import corner_moravec + corner_peaks, corner_fast, structure_tensor, + structure_tensor_eigvals, hessian_matrix, + hessian_matrix_eigvals) +from .corner_cy import corner_moravec, corner_orientations from .template import match_template +from .brief import BRIEF +from .censure import CENSURE +from .orb import ORB +from .match import match_descriptors +from .util import plot_matches __all__ = ['daisy', @@ -15,6 +22,10 @@ __all__ = ['daisy', 'greycoprops', 'local_binary_pattern', 'peak_local_max', + 'structure_tensor', + 'structure_tensor_eigvals', + 'hessian_matrix', + 'hessian_matrix_eigvals', 'corner_kitchen_rosenfeld', 'corner_harris', 'corner_shi_tomasi', @@ -22,4 +33,11 @@ __all__ = ['daisy', 'corner_subpix', 'corner_peaks', 'corner_moravec', - 'match_template'] + 'corner_fast', + 'corner_orientations', + 'match_template', + 'BRIEF', + 'CENSURE', + 'ORB', + 'match_descriptors', + 'plot_matches'] diff --git a/skimage/feature/_brief.py b/skimage/feature/_brief.py deleted file mode 100644 index 8a8d78db..00000000 --- a/skimage/feature/_brief.py +++ /dev/null @@ -1,229 +0,0 @@ -import numpy as np -from scipy.ndimage.filters import gaussian_filter - -from ..util import img_as_float -from .util import _mask_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): - """**Experimental function**. - - Extract BRIEF Descriptor about given keypoints for a given image. - - Parameters - ---------- - image : 2D ndarray - Input image. - keypoints : (P, 2) ndarray - 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. - 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 - Location i.e. (row, col) of keypoints after removing out those that - are near border. - - 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 - -------- - >> from skimage.feature import corner_peaks, corner_harris, \\ - .. 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 = keypoints[_mask_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): - """**Experimental function**. - - 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.py b/skimage/feature/brief.py new file mode 100644 index 00000000..d1626f17 --- /dev/null +++ b/skimage/feature/brief.py @@ -0,0 +1,181 @@ +import numpy as np +from scipy.ndimage.filters import gaussian_filter + +from .util import (DescriptorExtractor, _mask_border_keypoints, + _prepare_grayscale_input_2D) + +from .brief_cy import _brief_loop + + +class BRIEF(DescriptorExtractor): + + """BRIEF binary descriptor extractor. + + BRIEF (Binary Robust Independent Elementary Features) is an efficient + feature point descriptor. It is highly discriminative even when using + relatively few bits and is computed using simple intensity difference + tests. + + For each keypoint, intensity comparisons are carried out for a specifically + distributed number N of pixel-pairs resulting in a binary descriptor of + length N. For binary descriptors the Hamming distance can be used for + feature matching, which leads to lower computational cost in comparison to + the L2 norm. + + Parameters + ---------- + descriptor_size : int, optional + Size of BRIEF descriptor for each keypoint. Sizes 128, 256 and 512 + recommended by the authors. Default is 256. + patch_size : int, optional + Length of the two dimensional square patch sampling region around + the keypoints. Default is 49. + mode : {'normal', 'uniform'}, optional + Probability distribution for sampling location of decision pixel-pairs + around keypoints. + sample_seed : int, optional + Seed for the random sampling of 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` must be the same for the images to be + matched while building the descriptors. + sigma : float, optional + Standard deviation of the Gaussian low-pass filter applied to the image + to alleviate noise sensitivity, which is strongly recommended to obtain + discriminative and good descriptors. + + Attributes + ---------- + descriptors : (Q, `descriptor_size`) array of dtype bool + 2D ndarray of binary descriptors of size `descriptor_size` for Q + keypoints after filtering out border keypoints with value at an + index ``(i, j)`` either being ``True`` or ``False`` representing + the outcome of the intensity comparison for i-th keypoint on j-th + decision pixel-pair. It is ``Q == np.sum(mask)``. + mask : (N, ) array of dtype bool + Mask indicating whether a keypoint has been filtered out + (``False``) or is described in the `descriptors` array (``True``). + + Examples + -------- + >>> from skimage.feature import (corner_harris, corner_peaks, BRIEF, + ... match_descriptors) + >>> import numpy as np + >>> 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) + >>> 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) + >>> keypoints1 = corner_peaks(corner_harris(square1), min_distance=1) + >>> keypoints2 = corner_peaks(corner_harris(square2), min_distance=1) + >>> extractor = BRIEF(patch_size=5) + >>> extractor.extract(square1, keypoints1) + >>> descriptors1 = extractor.descriptors + >>> extractor.extract(square2, keypoints2) + >>> descriptors2 = extractor.descriptors + >>> matches = match_descriptors(descriptors1, descriptors2) + >>> matches + array([[0, 0], + [1, 1], + [2, 2], + [3, 3]]) + >>> keypoints1[matches[:, 0]] + array([[2, 2], + [2, 5], + [5, 2], + [5, 5]]) + >>> keypoints2[matches[:, 1]] + array([[2, 2], + [2, 6], + [6, 2], + [6, 6]]) + + """ + + def __init__(self, descriptor_size=256, patch_size=49, + mode='normal', sigma=1, sample_seed=1): + + mode = mode.lower() + if mode not in ('normal', 'uniform'): + raise ValueError("`mode` must be 'normal' or 'uniform'.") + + self.descriptor_size = descriptor_size + self.patch_size = patch_size + self.mode = mode + self.sigma = sigma + self.sample_seed = sample_seed + + self.descriptors = None + self.mask = None + + def extract(self, image, keypoints): + """Extract BRIEF binary descriptors for given keypoints in image. + + Parameters + ---------- + image : 2D array + Input image. + keypoints : (N, 2) array + Keypoint coordinates as ``(row, col)``. + + """ + + np.random.seed(self.sample_seed) + + image = _prepare_grayscale_input_2D(image) + + # Gaussian low-pass filtering to alleviate noise sensitivity + image = np.ascontiguousarray(gaussian_filter(image, self.sigma)) + + # Sampling pairs of decision pixels in patch_size x patch_size window + desc_size = self.descriptor_size + patch_size = self.patch_size + if self.mode == 'normal': + samples = (patch_size / 5.0) * np.random.randn(desc_size * 8) + samples = np.array(samples, dtype=np.int32) + samples = samples[(samples < (patch_size // 2)) + & (samples > - (patch_size - 2) // 2)] + + pos1 = samples[:desc_size * 2].reshape(desc_size, 2) + pos2 = samples[desc_size * 2:desc_size * 4].reshape(desc_size, 2) + elif self.mode == 'uniform': + samples = np.random.randint(-(patch_size - 2) // 2, + (patch_size // 2) + 1, + (desc_size * 2, 2)) + samples = np.array(samples, dtype=np.int32) + pos1, pos2 = np.split(samples, 2) + + pos1 = np.ascontiguousarray(pos1) + pos2 = np.ascontiguousarray(pos2) + + # Removing keypoints that are within (patch_size / 2) distance from the + # image border + self.mask = _mask_border_keypoints(image.shape, keypoints, + patch_size // 2) + + keypoints = np.array(keypoints[self.mask, :], dtype=np.intp, + order='C', copy=False) + + self.descriptors = np.zeros((keypoints.shape[0], desc_size), + dtype=bool, order='C') + + _brief_loop(image, self.descriptors.view(np.uint8), keypoints, + pos1, pos2) diff --git a/skimage/feature/_brief_cy.pyx b/skimage/feature/brief_cy.pyx similarity index 89% rename from skimage/feature/_brief_cy.pyx rename to skimage/feature/brief_cy.pyx index c53d85fc..8cd1afa7 100644 --- a/skimage/feature/_brief_cy.pyx +++ b/skimage/feature/brief_cy.pyx @@ -6,7 +6,7 @@ cimport numpy as cnp -def _brief_loop(double[:, ::1] image, char[:, ::1] descriptors, +def _brief_loop(double[:, ::1] image, unsigned char[:, ::1] descriptors, Py_ssize_t[:, ::1] keypoints, int[:, ::1] pos0, int[:, ::1] pos1): diff --git a/skimage/feature/censure.py b/skimage/feature/censure.py index 4bb7fdda..24d7e62d 100644 --- a/skimage/feature/censure.py +++ b/skimage/feature/censure.py @@ -1,9 +1,10 @@ import numpy as np from scipy.ndimage.filters import maximum_filter, minimum_filter, convolve +from skimage.feature.util import FeatureDetector, _prepare_grayscale_input_2D + from skimage.transform import integral_image -from skimage.feature.corner import _compute_auto_correlation -from skimage.util import img_as_float +from skimage.feature import structure_tensor from skimage.morphology import octagon, star from skimage.feature.util import _mask_border_keypoints @@ -65,19 +66,19 @@ def _filter_image(image, min_scale, max_scale, mode): 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)) + _octagon_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)) + response[:, :, i] = convolve(image, _star_kernel(m, n)) return response -def _octagon_filter_kernel(mo, no, mi, ni): +def _octagon_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) @@ -91,7 +92,7 @@ def _octagon_filter_kernel(mo, no, mi, ni): return bfilter -def _star_filter_kernel(m, n): +def _star_kernel(m, n): c = m + m // 2 - n - n // 2 outer_star = star(m) inner_star = np.zeros_like(outer_star) @@ -104,29 +105,25 @@ def _star_filter_kernel(m, n): 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 + Axx, Axy, Ayy = structure_tensor(image, sigma) + feature_mask[(Axx + Ayy) ** 2 + > line_threshold * (Axx * Ayy - Axy ** 2)] = False -def keypoints_censure(image, min_scale=1, max_scale=7, mode='DoB', - non_max_threshold=0.15, line_threshold=10): - """**Experimental function**. - Extracts CenSurE keypoints along with the corresponding scale using - either Difference of Boxes, Octagon or STAR bi-level filter. +class CENSURE(FeatureDetector): - Parameters - ---------- - image : 2D ndarray - Input image. - min_scale : int + """CENSURE keypoint detector. + + min_scale : int, optional Minimum scale to extract keypoints from. - max_scale : int + max_scale : int, optional 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'} + from the scales in the range [min_scale + 1, max_scale - 1]. The filter + sizes for different scales is such that the two adjacent scales + comprise of an octave. + mode : {'DoB', 'Octagon', 'STAR'}, optional 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 @@ -135,24 +132,24 @@ def keypoints_censure(image, min_scale=1, max_scale=7, mode='DoB', 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 + non_max_threshold : float, optional Threshold value used to suppress maximas and minimas with a weak magnitude response obtained after Non-Maximal Suppression. - line_threshold : float + line_threshold : float, optional Threshold for rejecting interest points which have ratio of principal curvatures greater than this value. - Returns - ------- + Attributes + ---------- 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. + Keypoint coordinates as ``(row, col)``. + scales : (N, ) array + Corresponding scales. References ---------- .. [1] Motilal Agrawal, Kurt Konolige and Morten Rufus Blas - "CenSurE: Center Surround Extremas for Realtime Feature + "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 @@ -161,74 +158,129 @@ def keypoints_censure(image, min_scale=1, max_scale=7, mode='DoB', Descriptors in the Context of Robot Navigation" http://www.jamris.org/01_2013/saveas.php?QUEST=JAMRIS_No01_2013_P_11-20.pdf + Examples + -------- + >>> from skimage.data import lena + >>> from skimage.color import rgb2gray + >>> from skimage.feature import CENSURE + >>> img = rgb2gray(lena()[100:300, 100:300]) + >>> censure = CENSURE() + >>> censure.detect(img) + >>> censure.keypoints + array([[ 71, 148], + [ 77, 186], + [ 78, 189], + [ 89, 174], + [127, 134], + [131, 133], + [134, 125], + [137, 125], + [149, 36], + [162, 165], + [168, 167], + [170, 5], + [171, 29], + [179, 20], + [194, 65]]) + >>> censure.scales + array([2, 4, 2, 3, 4, 2, 2, 3, 4, 6, 3, 2, 3, 4, 2]) + """ - # (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. + def __init__(self, min_scale=1, max_scale=7, mode='DoB', + non_max_threshold=0.15, line_threshold=10): - 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'.") - 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.') - 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.') + self.min_scale = min_scale + self.max_scale = max_scale + self.mode = mode + self.non_max_threshold = non_max_threshold + self.line_threshold = line_threshold - image = img_as_float(image) - image = np.ascontiguousarray(image) + self.keypoints = None + self.scales = None - # Generating all the scales - filter_response = _filter_image(image, min_scale, max_scale, mode) + def detect(self, image): + """Detect CENSURE keypoints along with the corresponding scale. - # 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 + Parameters + ---------- + image : 2D ndarray + Input image. - 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) + # (1) First we generate the required scales on the input grayscale + # image using a bi-level filter and stack them up in `filter_response`. - rows, cols, scales = np.nonzero(feature_mask[..., 1:max_scale - min_scale]) - keypoints = np.column_stack([rows, cols]) - scales = scales + min_scale + 1 + # (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. - if mode == 'dob': - return keypoints, scales + num_scales = self.max_scale - self.min_scale - cumulative_mask = np.zeros(keypoints.shape[0], dtype=np.bool) + image = np.ascontiguousarray(_prepare_grayscale_input_2D(image)) - 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) + # Generating all the scales + filter_response = _filter_image(image, self.min_scale, self.max_scale, + self.mode) - return keypoints[cumulative_mask], scales[cumulative_mask] + # Suppressing points that are neither minima or maxima in their + # 3 x 3 x 3 neighborhood 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 < self.non_max_threshold] = False + + for i in range(1, num_scales): + # 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 + (self.min_scale + i - 1) / 3.0), + self.line_threshold) + + rows, cols, scales = np.nonzero(feature_mask[..., 1:num_scales]) + keypoints = np.column_stack([rows, cols]) + scales = scales + self.min_scale + 1 + + if self.mode == 'dob': + self.keypoints = keypoints + self.scales = scales + return + + cumulative_mask = np.zeros(keypoints.shape[0], dtype=np.bool) + + if self.mode == 'octagon': + for i in range(self.min_scale + 1, self.max_scale): + c = (OCTAGON_OUTER_SHAPE[i - 1][0] - 1) // 2 \ + + OCTAGON_OUTER_SHAPE[i - 1][1] + cumulative_mask |= ( + _mask_border_keypoints(image.shape, keypoints, c) + & (scales == i)) + elif self.mode == 'star': + for i in range(self.min_scale + 1, self.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.shape, keypoints, c) + & (scales == i)) + + self.keypoints = keypoints[cumulative_mask] + self.scales = scales[cumulative_mask] diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index a025a7e4..48f74c27 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -1,18 +1,26 @@ import numpy as np from scipy import ndimage from scipy import stats + from skimage.color import rgb2grey from skimage.util import img_as_float, pad from skimage.feature import peak_local_max +from skimage.feature.util import _prepare_grayscale_input_2D +from skimage.feature.corner_cy import _corner_fast -def _compute_derivatives(image): +def _compute_derivatives(image, mode='constant', cval=0): """Compute derivatives in x and y direction using the Sobel operator. Parameters ---------- image : ndarray Input image. + mode : {'constant', 'reflect', 'wrap', 'nearest', 'mirror'}, optional + How to handle values outside the image borders. + cval : float, optional + Used in conjunction with mode 'constant', the value outside + the image boundaries. Returns ------- @@ -23,14 +31,82 @@ def _compute_derivatives(image): """ - imy = ndimage.sobel(image, axis=0, mode='constant', cval=0) - imx = ndimage.sobel(image, axis=1, mode='constant', cval=0) + imy = ndimage.sobel(image, axis=0, mode=mode, cval=cval) + imx = ndimage.sobel(image, axis=1, mode=mode, cval=cval) return imx, imy -def _compute_auto_correlation(image, sigma): - """Compute auto-correlation matrix using sum of squared differences. +def structure_tensor(image, sigma=1, mode='constant', cval=0): + """Compute structure tensor using sum of squared differences. + + The structure tensor A is defined as:: + + A = [Axx Axy] + [Axy Ayy] + + which is approximated by the weighted sum of squared differences in a local + window around each pixel in the image. + + Parameters + ---------- + image : ndarray + Input image. + sigma : float + Standard deviation used for the Gaussian kernel, which is used as a + weighting function for the local summation of squared differences. + mode : {'constant', 'reflect', 'wrap', 'nearest', 'mirror'}, optional + How to handle values outside the image borders. + cval : float, optional + Used in conjunction with mode 'constant', the value outside + the image boundaries. + + Returns + ------- + Axx : ndarray + Element of the structure tensor for each pixel in the input image. + Axy : ndarray + Element of the structure tensor for each pixel in the input image. + Ayy : ndarray + Element of the structure tensor for each pixel in the input image. + + Examples + -------- + >>> from skimage.feature import structure_tensor + >>> square = np.zeros((5, 5)) + >>> square[2, 2] = 1 + >>> Axx, Axy, Ayy = structure_tensor(square, sigma=0.1) + >>> Axx + array([[ 0., 0., 0., 0., 0.], + [ 0., 1., 0., 1., 0.], + [ 0., 4., 0., 4., 0.], + [ 0., 1., 0., 1., 0.], + [ 0., 0., 0., 0., 0.]]) + + """ + + image = _prepare_grayscale_input_2D(image) + + imx, imy = _compute_derivatives(image, mode=mode, cval=cval) + + # structure tensore + Axx = ndimage.gaussian_filter(imx * imx, sigma, mode=mode, cval=cval) + Axy = ndimage.gaussian_filter(imx * imy, sigma, mode=mode, cval=cval) + Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode=mode, cval=cval) + + return Axx, Axy, Ayy + + +def hessian_matrix(image, sigma=1, mode='constant', cval=0): + """Compute Hessian matrix. + + The Hessian matrix is defined as:: + + H = [Hxx Hxy] + [Hxy Hyy] + + which is computed by convolving the image with the second derivatives + of the Gaussian kernel in the respective x- and y-directions. Parameters ---------- @@ -39,32 +115,142 @@ def _compute_auto_correlation(image, sigma): sigma : float Standard deviation used for the Gaussian kernel, which is used as weighting function for the auto-correlation matrix. + mode : {'constant', 'reflect', 'wrap', 'nearest', 'mirror'}, optional + How to handle values outside the image borders. + cval : float, optional + Used in conjunction with mode 'constant', the value outside + the image boundaries. Returns ------- - Axx : ndarray - Element of the auto-correlation matrix for each pixel in input image. - Axy : ndarray - Element of the auto-correlation matrix for each pixel in input image. - Ayy : ndarray - Element of the auto-correlation matrix for each pixel in input image. + Hxx : ndarray + Element of the Hessian matrix for each pixel in the input image. + Hxy : ndarray + Element of the Hessian matrix for each pixel in the input image. + Hyy : ndarray + Element of the Hessian matrix for each pixel in the input image. + + Examples + -------- + >>> from skimage.feature import hessian_matrix, hessian_matrix_eigvals + >>> square = np.zeros((5, 5)) + >>> square[2, 2] = 1 + >>> Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) + >>> Hxx + array([[ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 1., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.]]) """ - if image.ndim == 3: - image = img_as_float(rgb2grey(image)) + image = _prepare_grayscale_input_2D(image) - imx, imy = _compute_derivatives(image) + # window extent to the left and right, which covers > 99% of the normal + # distribution + window_ext = max(1, np.ceil(3 * sigma)) - # structure tensore - Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0) - Axy = ndimage.gaussian_filter(imx * imy, sigma, mode='constant', cval=0) - Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode='constant', cval=0) + ky, kx = np.mgrid[-window_ext:window_ext + 1, -window_ext:window_ext + 1] - return Axx, Axy, Ayy + # second derivative Gaussian kernels + gaussian_exp = np.exp(-(kx ** 2 + ky ** 2) / (2 * sigma ** 2)) + kernel_xx = 1 / (2 * np.pi * sigma ** 4) * (kx ** 2 / sigma ** 2 - 1) + kernel_xx *= gaussian_exp + kernel_xx /= kernel_xx.sum() + kernel_xy = 1 / (2 * np.pi * sigma ** 6) * (kx * ky) + kernel_xy *= gaussian_exp + kernel_xy /= kernel_xx.sum() + kernel_yy = kernel_xx.transpose() + + Hxx = ndimage.convolve(image, kernel_xx, mode=mode, cval=cval) + Hxy = ndimage.convolve(image, kernel_xy, mode=mode, cval=cval) + Hyy = ndimage.convolve(image, kernel_yy, mode=mode, cval=cval) + + return Hxx, Hxy, Hyy -def corner_kitchen_rosenfeld(image): +def _image_orthogonal_matrix22_eigvals(M00, M01, M11): + l1 = (M00 + M11) / 2 + np.sqrt(4 * M01 ** 2 + (M00 - M11) ** 2) / 2 + l2 = (M00 + M11) / 2 - np.sqrt(4 * M01 ** 2 + (M00 - M11) ** 2) / 2 + return l1, l2 + + +def structure_tensor_eigvals(Axx, Axy, Ayy): + """Compute Eigen values of structure tensor. + + Parameters + ---------- + Axx : ndarray + Element of the structure tensor for each pixel in the input image. + Axy : ndarray + Element of the structure tensor for each pixel in the input image. + Ayy : ndarray + Element of the structure tensor for each pixel in the input image. + + Returns + ------- + l1 : ndarray + Larger eigen value for each input matrix. + l2 : ndarray + Smaller eigen value for each input matrix. + + Examples + -------- + >>> from skimage.feature import structure_tensor, structure_tensor_eigvals + >>> square = np.zeros((5, 5)) + >>> square[2, 2] = 1 + >>> Axx, Axy, Ayy = structure_tensor(square, sigma=0.1) + >>> structure_tensor_eigvals(Axx, Axy, Ayy)[0] + array([[ 0., 0., 0., 0., 0.], + [ 0., 2., 4., 2., 0.], + [ 0., 4., 0., 4., 0.], + [ 0., 2., 4., 2., 0.], + [ 0., 0., 0., 0., 0.]]) + + """ + + return _image_orthogonal_matrix22_eigvals(Axx, Axy, Ayy) + + +def hessian_matrix_eigvals(Hxx, Hxy, Hyy): + """Compute Eigen values of Hessian matrix. + + Parameters + ---------- + Hxx : ndarray + Element of the Hessian matrix for each pixel in the input image. + Hxy : ndarray + Element of the Hessian matrix for each pixel in the input image. + Hyy : ndarray + Element of the Hessian matrix for each pixel in the input image. + + Returns + ------- + l1 : ndarray + Larger eigen value for each input matrix. + l2 : ndarray + Smaller eigen value for each input matrix. + + Examples + -------- + >>> from skimage.feature import hessian_matrix, hessian_matrix_eigvals + >>> square = np.zeros((5, 5)) + >>> square[2, 2] = 1 + >>> Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) + >>> hessian_matrix_eigvals(Hxx, Hxy, Hyy)[0] + array([[ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 1., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.]]) + + """ + + return _image_orthogonal_matrix22_eigvals(Hyy, Hxy, Hyy) + + +def corner_kitchen_rosenfeld(image, mode='constant', cval=0): """Compute Kitchen and Rosenfeld corner measure response image. The corner measure is calculated as follows:: @@ -79,6 +265,11 @@ def corner_kitchen_rosenfeld(image): ---------- image : ndarray Input image. + mode : {'constant', 'reflect', 'wrap', 'nearest', 'mirror'}, optional + How to handle values outside the image borders. + cval : float, optional + Used in conjunction with mode 'constant', the value outside + the image boundaries. Returns ------- @@ -87,9 +278,9 @@ def corner_kitchen_rosenfeld(image): """ - imx, imy = _compute_derivatives(image) - imxx, imxy = _compute_derivatives(imx) - imyx, imyy = _compute_derivatives(imy) + imx, imy = _compute_derivatives(image, mode=mode, cval=cval) + imxx, imxy = _compute_derivatives(imx, mode=mode, cval=cval) + imyx, imyy = _compute_derivatives(imy, mode=mode, cval=cval) numerator = (imxx * imy**2 + imyy * imx**2 - 2 * imxy * imx * imy) denominator = (imx**2 + imy**2) @@ -147,9 +338,9 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): Examples -------- >>> from skimage.feature import corner_harris, corner_peaks - >>> square = np.zeros([10, 10], dtype=int) + >>> square = np.zeros([10, 10]) >>> square[2:8, 2:8] = 1 - >>> square + >>> square.astype(int) array([[0, 0, 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, 1, 0, 0], @@ -168,7 +359,7 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """ - Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) + Axx, Axy, Ayy = structure_tensor(image, sigma) # determinant detA = Axx * Ayy - Axy**2 @@ -217,9 +408,9 @@ def corner_shi_tomasi(image, sigma=1): Examples -------- >>> from skimage.feature import corner_shi_tomasi, corner_peaks - >>> square = np.zeros([10, 10], dtype=int) + >>> square = np.zeros([10, 10]) >>> square[2:8, 2:8] = 1 - >>> square + >>> square.astype(int) array([[0, 0, 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, 1, 0, 0], @@ -238,7 +429,7 @@ def corner_shi_tomasi(image, sigma=1): """ - Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) + Axx, Axy, Ayy = structure_tensor(image, sigma) # minimum eigenvalue of A response = ((Axx + Ayy) - np.sqrt((Axx - Ayy)**2 + 4 * Axy**2)) / 2 @@ -308,7 +499,7 @@ def corner_foerstner(image, sigma=1): """ - Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) + Axx, Axy, Ayy = structure_tensor(image, sigma) # determinant detA = Axx * Ayy - Axy**2 @@ -326,6 +517,69 @@ def corner_foerstner(image, sigma=1): return w, q +def corner_fast(image, n=12, threshold=0.15): + """Extract FAST corners for a given image. + + Parameters + ---------- + image : 2D ndarray + Input image. + n : int + Minimum number of consecutive pixels out of 16 pixels on the circle + that should all be either brighter or darker w.r.t testpixel. + A point c on the circle is darker w.r.t test pixel p if + `Ic < Ip - threshold` and brighter if `Ic > Ip + threshold`. Also + stands for the n in `FAST-n` corner detector. + threshold : float + Threshold used in deciding whether the pixels on the circle are + brighter, darker or similar w.r.t. the test pixel. Decrease the + threshold when more corners are desired and vice-versa. + + Returns + ------- + response : ndarray + FAST corner response image. + + References + ---------- + .. [1] Edward Rosten and Tom Drummond + "Machine Learning for high-speed corner detection", + http://www.edwardrosten.com/work/rosten_2006_machine.pdf + .. [2] Wikipedia, "Features from accelerated segment test", + https://en.wikipedia.org/wiki/Features_from_accelerated_segment_test + + Examples + -------- + >>> from skimage.feature import corner_fast, corner_peaks + >>> square = np.zeros((12, 12)) + >>> square[3:9, 3:9] = 1 + >>> square.astype(int) + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 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, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 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, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) + >>> corner_peaks(corner_fast(square, 9), min_distance=1) + array([[3, 3], + [3, 8], + [8, 3], + [8, 8]]) + + """ + image = _prepare_grayscale_input_2D(image) + + image = np.ascontiguousarray(image) + response = _corner_fast(image, n, threshold) + return response + + def corner_subpix(image, corners, window_size=11, alpha=0.99): """Determine subpixel position of corners. @@ -354,10 +608,10 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): Examples -------- >>> from skimage.feature import corner_harris, corner_peaks, corner_subpix - >>> img = np.zeros((10, 10), dtype=int) + >>> img = np.zeros((10, 10)) >>> img[:5, :5] = 1 >>> img[5:, 5:] = 1 - >>> img + >>> img.astype(int) array([[1, 1, 1, 1, 1, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 0, 0, 0, 0, 0], @@ -408,7 +662,7 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): maxx = x0 + wext + 2 window = image[miny:maxy, minx:maxx] - winx, winy = _compute_derivatives(window) + winx, winy = _compute_derivatives(window, mode='constant', cval=0) # compute gradient suares and remove border winx_winx = (winx * winx)[1:-1, 1:-1] diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx index c459ee92..86ad3c43 100644 --- a/skimage/feature/corner_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -5,9 +5,12 @@ import numpy as np cimport numpy as cnp from libc.float cimport DBL_MAX +from libc.math cimport atan2 +from skimage.util import img_as_float, pad from skimage.color import rgb2grey -from skimage.util import img_as_float + +from .util import _prepare_grayscale_input_2D def corner_moravec(image, Py_ssize_t window_size=1): @@ -30,30 +33,30 @@ def corner_moravec(image, Py_ssize_t window_size=1): References ---------- - ..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/moravec.htm - ..[2] http://en.wikipedia.org/wiki/Corner_detection + .. [1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/moravec.htm + .. [2] http://en.wikipedia.org/wiki/Corner_detection Examples -------- - >>> from skimage.feature import corner_moravec, peak_local_max + >>> from skimage.feature import corner_moravec >>> square = np.zeros([7, 7]) >>> square[3, 3] = 1 - >>> square - array([[ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 1., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.]]) - >>> corner_moravec(square) - array([[ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 1., 1., 1., 0., 0.], - [ 0., 0., 1., 2., 1., 0., 0.], - [ 0., 0., 1., 1., 1., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0.]]) + >>> square.astype(int) + array([[0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0]]) + >>> corner_moravec(square).astype(int) + array([[0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 1, 0, 0], + [0, 0, 1, 2, 1, 0, 0], + [0, 0, 1, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0]]) """ cdef Py_ssize_t rows = image.shape[0] @@ -80,3 +83,193 @@ def corner_moravec(image, Py_ssize_t window_size=1): out[r, c] = min_msum return np.asarray(out) + + +cdef inline double _corner_fast_response(double curr_pixel, + double* circle_intensities, + char* bins, char state, char n): + cdef char consecutive_count = 0 + cdef double curr_response + cdef Py_ssize_t l, m + for l in range(15 + n): + if bins[l % 16] == state: + consecutive_count += 1 + if consecutive_count == n: + curr_response = 0 + for m in range(16): + curr_response += abs(circle_intensities[m] - curr_pixel) + return curr_response + else: + consecutive_count = 0 + return 0 + + +def _corner_fast(double[:, ::1] image, char n, double threshold): + + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + + cdef Py_ssize_t i, j, k + + cdef char speed_sum_b, speed_sum_d + cdef double curr_pixel + cdef double lower_threshold, upper_threshold + cdef double[:, ::1] corner_response = np.zeros((rows, cols), + dtype=np.double) + + cdef char *rp = [0, 1, 2, 3, 3, 3, 2, 1, 0, -1, -2, -3, -3, -3, -2, -1] + cdef char *cp = [3, 3, 2, 1, 0, -1, -2, -3, -3, -3, -2, -1, 0, 1, 2, 3] + cdef char bins[16] + cdef double circle_intensities[16] + + cdef double curr_response + + for i in range(3, rows - 3): + for j in range(3, cols - 3): + + curr_pixel = image[i, j] + lower_threshold = curr_pixel - threshold + upper_threshold = curr_pixel + threshold + + for k in range(16): + circle_intensities[k] = image[i + rp[k], j + cp[k]] + if circle_intensities[k] > upper_threshold: + # Brighter pixel + bins[k] = 'b' + elif circle_intensities[k] < lower_threshold: + # Darker pixel + bins[k] = 'd' + else: + # Similar pixel + bins[k] = 's' + + # High speed test for n >= 12 + if n >= 12: + speed_sum_b = 0 + speed_sum_d = 0 + for k in range(0, 16, 4): + if bins[k] == 'b': + speed_sum_b += 1 + elif bins[k] == 'd': + speed_sum_d += 1 + if speed_sum_d < 3 and speed_sum_b < 3: + continue + + # Test for bright pixels + curr_response = \ + _corner_fast_response(curr_pixel, circle_intensities, + bins, 'b', n) + + # Test for dark pixels + if curr_response == 0: + curr_response = \ + _corner_fast_response(curr_pixel, circle_intensities, + bins, 'd', n) + + corner_response[i, j] = curr_response + + return np.asarray(corner_response) + + +def corner_orientations(image, Py_ssize_t[:, :] corners, mask): + """Compute the orientation of corners. + + The orientation of corners is computed using the first order central moment + i.e. the center of mass approach. The corner orientation is the angle of + the vector from the corner coordinate to the intensity centroid in the + local neighborhood around the corner calculated using first order central + moment. + + Parameters + ---------- + image : 2D array + Input grayscale image. + corners : (N, 2) array + Corner coordinates as ``(row, col)``. + mask : 2D array + Mask defining the local neighborhood of the corner used for the + calculation of the central moment. + + Returns + ------- + orientations : (N, 1) array + Orientations of corners in the range [-pi, pi]. + + References + ---------- + .. [1] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski + "ORB : An efficient alternative to SIFT and SURF" + http://www.vision.cs.chubu.ac.jp/CV-R/pdf/Rublee_iccv2011.pdf + .. [2] Paul L. Rosin, "Measuring Corner Properties" + http://users.cs.cf.ac.uk/Paul.Rosin/corner2.pdf + + Examples + -------- + >>> from skimage.morphology import octagon + >>> from skimage.feature import (corner_fast, corner_peaks, + ... corner_orientations) + >>> square = np.zeros((12, 12)) + >>> square[3:9, 3:9] = 1 + >>> square.astype(int) + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 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, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 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, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) + >>> corners = corner_peaks(corner_fast(square, 9), min_distance=1) + >>> corners + array([[3, 3], + [3, 8], + [8, 3], + [8, 8]]) + >>> orientations = corner_orientations(square, corners, octagon(3, 2)) + >>> np.rad2deg(orientations) + array([ 45., 135., -45., -135.]) + + """ + + image = _prepare_grayscale_input_2D(image) + + if mask.shape[0] % 2 != 1 or mask.shape[1] % 2 != 1: + raise ValueError("Size of mask must be uneven.") + + cdef unsigned char[:, ::1] cmask = np.ascontiguousarray(mask != 0, + dtype=np.uint8) + + cdef Py_ssize_t i, r, c, r0, c0 + cdef Py_ssize_t mrows = mask.shape[0] + cdef Py_ssize_t mcols = mask.shape[1] + cdef Py_ssize_t mrows2 = (mrows - 1) / 2 + cdef Py_ssize_t mcols2 = (mcols - 1) / 2 + cdef double[:, :] cimage = pad(image, (mrows2, mcols2), mode='constant', + constant_values=0) + cdef double[:] orientations = np.zeros(corners.shape[0], dtype=np.double) + cdef double curr_pixel + cdef double m01, m10, m01_tmp + + for i in range(corners.shape[0]): + r0 = corners[i, 0] + c0 = corners[i, 1] + + m01 = 0 + m10 = 0 + + for r in range(mrows): + m01_tmp = 0 + for c in range(mcols): + if cmask[r, c]: + curr_pixel = cimage[r0 + r, c0 + c] + m10 += curr_pixel * (c - mcols2) + m01_tmp += curr_pixel + m01 += m01_tmp * (r - mrows2) + + orientations[i] = atan2(m01, m10) + + return np.asarray(orientations) diff --git a/skimage/feature/match.py b/skimage/feature/match.py new file mode 100644 index 00000000..ccd8c035 --- /dev/null +++ b/skimage/feature/match.py @@ -0,0 +1,65 @@ +import numpy as np +from scipy.spatial.distance import cdist + + +def match_descriptors(descriptors1, descriptors2, metric=None, p=2, + threshold=0, cross_check=True): + """Brute-force matching of descriptors. + + For each descriptor in the first set this matcher finds the closest + descriptor in the second set (and vice-versa in the case of enabled + cross-checking). + + Parameters + ---------- + descriptors1 : (M, P) array + Binary descriptors of size P about M keypoints in the first image. + descriptors2 : (N, P) array + Binary descriptors of size P about N keypoints in the second image. + metric : {'euclidean', 'cityblock', 'minkowski', 'hamming', ...} + The metric to compute the distance between two descriptors. See + `scipy.spatial.distance.cdist` for all possible types. The hamming + distance should be used for binary descriptors. By default the L2-norm + is used for all descriptors of dtype float or double and the Hamming + distance is used for binary descriptors automatically. + p : int + The p-norm to apply for ``metric='minkowski'``. + threshold : float + Maximum allowed distance between descriptors of two keypoints + in separate images to be regarded as a match. + cross_check : bool + If True, the matched keypoints are returned after cross checking i.e. a + matched pair (keypoint1, keypoint2) is returned if keypoint2 is the + best match for keypoint1 in second image and keypoint1 is the best + match for keypoint2 in first image. + + Returns + ------- + matches : (Q, 2) array + Indices of corresponding matches in first and second set of + descriptors, where ``matches[:, 0]`` denote the indices in the first + and ``matches[:, 1]`` the indices in the second set of descriptors. + + """ + + if descriptors1.shape[1] != descriptors2.shape[1]: + raise ValueError("Descriptor length must equal.") + + if metric is None: + if np.issubdtype(descriptors1.dtype, np.bool): + metric = 'hamming' + else: + metric = 'euclidean' + + distances = cdist(descriptors1, descriptors2, metric=metric, p=p) + + indices1 = np.arange(descriptors1.shape[0]) + indices2 = np.argmin(distances, axis=1) + + if cross_check: + matches1 = np.argmin(distances, axis=0) + mask = indices1 == matches1[indices2] + indices1 = indices1[mask] + indices2 = indices2[mask] + + return np.column_stack((indices1, indices2)) diff --git a/skimage/feature/orb.py b/skimage/feature/orb.py new file mode 100644 index 00000000..2ddcf4f3 --- /dev/null +++ b/skimage/feature/orb.py @@ -0,0 +1,336 @@ +import numpy as np + +from skimage.feature.util import (FeatureDetector, DescriptorExtractor, + _mask_border_keypoints, + _prepare_grayscale_input_2D) + +from skimage.feature import (corner_fast, corner_orientations, corner_peaks, + corner_harris) +from skimage.transform import pyramid_gaussian + +from .orb_cy import _orb_loop + + +OFAST_MASK = np.zeros((31, 31)) +OFAST_UMAX = [15, 15, 15, 15, 14, 14, 14, 13, 13, 12, 11, 10, 9, 8, 6, 3] +for i in range(-15, 16): + for j in range(-OFAST_UMAX[abs(i)], OFAST_UMAX[abs(i)] + 1): + OFAST_MASK[15 + j, 15 + i] = 1 + + +class ORB(FeatureDetector, DescriptorExtractor): + + """Oriented FAST and rotated BRIEF feature detector and binary descriptor + extractor. + + Parameters + ---------- + n_keypoints : int, optional + Number of keypoints to be returned. The function will return the best + `n_keypoints` according to the Harris corner response if more than + `n_keypoints` are detected. If not, then all the detected keypoints + are returned. + fast_n : int, optional + The `n` parameter in `skimage.feature.corner_fast`. Minimum number of + consecutive pixels out of 16 pixels on the circle that should all be + either brighter or darker w.r.t test-pixel. A point c on the circle is + darker w.r.t test pixel p if ``Ic < Ip - threshold`` and brighter if + ``Ic > Ip + threshold``. Also stands for the n in ``FAST-n`` corner + detector. + fast_threshold : float, optional + The ``threshold`` parameter in ``feature.corner_fast``. Threshold used + to decide whether the pixels on the circle are brighter, darker or + similar w.r.t. the test pixel. Decrease the threshold when more + corners are desired and vice-versa. + harris_k : float, optional + The `k` parameter in `skimage.feature.corner_harris`. Sensitivity + factor to separate corners from edges, typically in range ``[0, 0.2]``. + Small values of `k` result in detection of sharp corners. + downscale : float, optional + Downscale factor for the image pyramid. Default value 1.2 is chosen so + that there are more dense scales which enable robust scale invariance + for a subsequent feature description. + n_scales : int, optional + Maximum number of scales from the bottom of the image pyramid to + extract the features from. + + Attributes + ---------- + keypoints : (N, 2) array + Keypoint coordinates as ``(row, col)``. + scales : (N, ) array + Corresponding scales. + orientations : (N, ) array + Corresponding orientations in radians. + responses : (N, ) array + Corresponding Harris corner responses. + descriptors : (Q, `descriptor_size`) array of dtype bool + 2D array of binary descriptors of size `descriptor_size` for Q + keypoints after filtering out border keypoints with value at an + index ``(i, j)`` either being ``True`` or ``False`` representing + the outcome of the intensity comparison for i-th keypoint on j-th + decision pixel-pair. It is ``Q == np.sum(mask)``. + + References + ---------- + .. [1] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski + "ORB: An efficient alternative to SIFT and SURF" + http://www.vision.cs.chubu.ac.jp/CV-R/pdf/Rublee_iccv2011.pdf + + Examples + -------- + >>> from skimage.feature import ORB, match_descriptors + >>> img1 = np.zeros((100, 100)) + >>> img2 = np.zeros_like(img1) + >>> np.random.seed(1) + >>> square = np.random.rand(20, 20) + >>> img1[40:60, 40:60] = square + >>> img2[53:73, 53:73] = square + >>> detector_extractor1 = ORB(n_keypoints=5) + >>> detector_extractor2 = ORB(n_keypoints=5) + >>> detector_extractor1.detect_and_extract(img1) + >>> detector_extractor2.detect_and_extract(img2) + >>> matches = match_descriptors(detector_extractor1.descriptors, + ... detector_extractor2.descriptors) + >>> matches + array([[0, 0], + [1, 1], + [2, 2], + [3, 3], + [4, 4]]) + >>> detector_extractor1.keypoints[matches[:, 0]] + array([[ 42., 40.], + [ 47., 58.], + [ 44., 40.], + [ 59., 42.], + [ 45., 44.]]) + >>> detector_extractor2.keypoints[matches[:, 1]] + array([[ 55., 53.], + [ 60., 71.], + [ 57., 53.], + [ 72., 55.], + [ 58., 57.]]) + + """ + + def __init__(self, downscale=1.2, n_scales=8, + n_keypoints=500, fast_n=9, fast_threshold=0.08, + harris_k=0.04): + self.downscale = downscale + self.n_scales = n_scales + self.n_keypoints = n_keypoints + self.fast_n = fast_n + self.fast_threshold = fast_threshold + self.harris_k = harris_k + + self.keypoints = None + self.scales = None + self.responses = None + self.orientations = None + self.descriptors = None + + def _build_pyramid(self, image): + image = _prepare_grayscale_input_2D(image) + return list(pyramid_gaussian(image, self.n_scales - 1, self.downscale)) + + def _detect_octave(self, octave_image): + # Extract keypoints for current octave + fast_response = corner_fast(octave_image, self.fast_n, + self.fast_threshold) + keypoints = corner_peaks(fast_response, min_distance=1) + + if len(keypoints) == 0: + return (np.zeros((0, 2), dtype=np.double), + np.zeros((0, ), dtype=np.double), + np.zeros((0, ), dtype=np.double)) + + mask = _mask_border_keypoints(octave_image.shape, keypoints, + distance=16) + keypoints = keypoints[mask] + + orientations = corner_orientations(octave_image, keypoints, + OFAST_MASK) + + harris_response = corner_harris(octave_image, method='k', + k=self.harris_k) + responses = harris_response[keypoints[:, 0], keypoints[:, 1]] + + return keypoints, orientations, responses + + def detect(self, image): + """Detect oriented FAST keypoints along with the corresponding scale. + + Parameters + ---------- + image : 2D array + Input image. + + """ + + pyramid = self._build_pyramid(image) + + keypoints_list = [] + orientations_list = [] + scales_list = [] + responses_list = [] + + for octave in range(len(pyramid)): + + octave_image = np.ascontiguousarray(pyramid[octave]) + + keypoints, orientations, responses = \ + self._detect_octave(octave_image) + + keypoints_list.append(keypoints * self.downscale ** octave) + orientations_list.append(orientations) + scales_list.append(self.downscale ** octave + * np.ones(keypoints.shape[0], dtype=np.intp)) + responses_list.append(responses) + + keypoints = np.vstack(keypoints_list) + orientations = np.hstack(orientations_list) + scales = np.hstack(scales_list) + responses = np.hstack(responses_list) + + if keypoints.shape[0] < self.n_keypoints: + self.keypoints = keypoints + self.scales = scales + self.orientations = orientations + self.responses = responses + else: + # Choose best n_keypoints according to Harris corner response + best_indices = responses.argsort()[::-1][:self.n_keypoints] + self.keypoints = keypoints[best_indices] + self.scales = scales[best_indices] + self.orientations = orientations[best_indices] + self.responses = responses[best_indices] + + def _extract_octave(self, octave_image, keypoints, orientations): + mask = _mask_border_keypoints(octave_image.shape, keypoints, + distance=20) + keypoints = np.array(keypoints[mask], dtype=np.intp, order='C', + copy=False) + orientations = np.array(orientations[mask], dtype=np.double, order='C', + copy=False) + + descriptors = _orb_loop(octave_image, keypoints, orientations) + + return descriptors, mask + + def extract(self, image, keypoints, scales, orientations): + """Extract rBRIEF binary descriptors for given keypoints in image. + + Note that the keypoints must be extracted using the same `downscale` + and `n_scales` parameters. Additionally, if you want to extract both + keypoints and descriptors you should use the faster + `detect_and_extract`. + + Parameters + ---------- + image : 2D array + Input image. + keypoints : (N, 2) array + Keypoint coordinates as ``(row, col)``. + scales : (N, ) array + Corresponding scales. + orientations : (N, ) array + Corresponding orientations in radians. + + """ + + pyramid = self._build_pyramid(image) + + descriptors_list = [] + mask_list = [] + + # Determine octaves from scales + octaves = (np.log(scales) / np.log(self.downscale)).astype(np.intp) + + for octave in range(len(pyramid)): + + # Mask for all keypoints in current octave + octave_mask = octaves == octave + + if np.sum(octave_mask) > 0: + + octave_image = np.ascontiguousarray(pyramid[octave]) + + octave_keypoints = keypoints[octave_mask] + octave_keypoints /= self.downscale ** octave + + octave_orientations = orientations[octave_mask] + + descriptors, mask = self._extract_octave(octave_image, + octave_keypoints, + octave_orientations) + + descriptors_list.append(descriptors) + mask_list.append(mask) + + self.descriptors = np.vstack(descriptors_list).view(np.bool) + self.mask_ = np.hstack(mask_list) + + def detect_and_extract(self, image): + """Detect oriented FAST keypoints and extract rBRIEF descriptors. + + Note that this is faster than first calling `detect` and then + `extract`. + + Parameters + ---------- + image : 2D array + Input image. + + """ + + pyramid = self._build_pyramid(image) + + keypoints_list = [] + responses_list = [] + scales_list = [] + orientations_list = [] + descriptors_list = [] + + for octave in range(len(pyramid)): + + octave_image = np.ascontiguousarray(pyramid[octave]) + + keypoints, orientations, responses = \ + self._detect_octave(octave_image) + + if len(keypoints) == 0: + keypoints_list.append(keypoints) + responses_list.append(responses) + descriptors_list.append(np.zeros((0, 256), dtype=np.bool)) + continue + + descriptors, mask = self._extract_octave(octave_image, keypoints, + orientations) + + keypoints_list.append(keypoints[mask] * self.downscale ** octave) + responses_list.append(responses[mask]) + orientations_list.append(orientations[mask]) + scales_list.append(self.downscale ** octave + * np.ones(keypoints.shape[0], dtype=np.intp)) + descriptors_list.append(descriptors) + + keypoints = np.vstack(keypoints_list) + responses = np.hstack(responses_list) + scales = np.hstack(scales_list) + orientations = np.hstack(orientations_list) + descriptors = np.vstack(descriptors_list).view(np.bool) + + if keypoints.shape[0] < self.n_keypoints: + self.keypoints = keypoints + self.scales = scales + self.orientations = orientations + self.responses = responses + self.descriptors = descriptors + else: + # Choose best n_keypoints according to Harris corner response + best_indices = responses.argsort()[::-1][:self.n_keypoints] + self.keypoints = keypoints[best_indices] + self.scales = scales[best_indices] + self.orientations = orientations[best_indices] + self.responses = responses[best_indices] + self.descriptors = descriptors[best_indices] diff --git a/skimage/feature/orb_cy.pyx b/skimage/feature/orb_cy.pyx new file mode 100644 index 00000000..b497c74d --- /dev/null +++ b/skimage/feature/orb_cy.pyx @@ -0,0 +1,54 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import os +import numpy as np + +from skimage import data_dir + +cimport numpy as cnp +from libc.math cimport sin, cos, round + +POS = np.loadtxt(os.path.join(data_dir, "orb_descriptor_positions.txt"), + dtype=np.int8) +POS0 = np.ascontiguousarray(POS[:, :2]) +POS1 = np.ascontiguousarray(POS[:, 2:]) + + +def _orb_loop(double[:, ::1] image, Py_ssize_t[:, ::1] keypoints, + double[:] orientations): + + cdef Py_ssize_t i, d, kr, kc, pr0, pr1, pc0, pc1, spr0, spc0, spr1, spc1 + cdef int[:, ::1] steered_pos0, steered_pos1 + cdef double angle + cdef char[:, ::1] descriptors = np.zeros((keypoints.shape[0], + POS.shape[0]), dtype=np.uint8) + cdef char[:, ::1] cpos0 = POS0 + cdef char[:, ::1] cpos1 = POS1 + + for i in range(descriptors.shape[0]): + + angle = orientations[i] + sin_a = sin(angle) + cos_a = cos(angle) + + kr = keypoints[i, 0] + kc = keypoints[i, 1] + + for j in range(descriptors.shape[1]): + pr0 = cpos0[j, 0] + pc0 = cpos0[j, 1] + pr1 = cpos1[j, 0] + pc1 = cpos1[j, 1] + + spr0 = round(sin_a * pr0 + cos_a * pc0) + spc0 = round(cos_a * pr0 - sin_a * pc0) + spr1 = round(sin_a * pr1 + cos_a * pc1) + spc1 = round(cos_a * pr1 - sin_a * pc1) + + if image[kr + spr0, kc + spc0] < image[kr + spr1, kc + spc1]: + descriptors[i, j] = True + + return np.asarray(descriptors) diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 915bc351..9d8b4fbb 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -14,14 +14,17 @@ def configuration(parent_package='', top_path=None): 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(['orb_cy.pyx'], working_path=base_path) + cython(['brief_cy.pyx'], working_path=base_path) cython(['_texture.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'], + config.add_extension('orb_cy', sources=['orb_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']) diff --git a/skimage/feature/tests/_test_brief.py b/skimage/feature/tests/_test_brief.py deleted file mode 100644 index 1d26cbbd..00000000 --- a/skimage/feature/tests/_test_brief.py +++ /dev/null @@ -1,83 +0,0 @@ -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.color import rgb2gray -from skimage.feature import (brief, match_keypoints_brief, corner_peaks, - corner_harris) - - -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_censure.py b/skimage/feature/tests/_test_censure.py deleted file mode 100644 index 4cd2ad68..00000000 --- a/skimage/feature/tests/_test_censure.py +++ /dev/null @@ -1,89 +0,0 @@ -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/tests/test_brief.py b/skimage/feature/tests/test_brief.py new file mode 100644 index 00000000..554301b2 --- /dev/null +++ b/skimage/feature/tests/test_brief.py @@ -0,0 +1,77 @@ +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.color import rgb2gray +from skimage.feature import BRIEF, corner_peaks, corner_harris + + +def test_color_image_unsupported_error(): + """Brief descriptors can be evaluated on gray-scale images only.""" + img = np.zeros((20, 20, 3)) + keypoints = np.asarray([[7, 5], [11, 13]]) + assert_raises(ValueError, BRIEF().extract, img, keypoints) + + +def test_normal_mode(): + """Verify the computed BRIEF descriptors with expected for normal mode.""" + img = rgb2gray(data.lena()) + + keypoints = corner_peaks(corner_harris(img), min_distance=5) + + extractor = BRIEF(descriptor_size=8, sigma=2) + + extractor.extract(img, keypoints[:8]) + + expected = np.array([[ True, False, True, False, True, True, False, False], + [False, False, False, False, True, False, False, False], + [ True, True, True, True, True, True, True, True], + [ True, False, True, True, False, True, False, True], + [False, True, True, True, True, True, True, True], + [ True, False, False, False, False, True, False, True], + [False, True, True, True, False, False, True, False], + [False, False, False, False, True, False, False, False]], dtype=bool) + + assert_array_equal(extractor.descriptors, expected) + + +def test_uniform_mode(): + """Verify the computed BRIEF descriptors with expected for uniform mode.""" + img = rgb2gray(data.lena()) + + keypoints = corner_peaks(corner_harris(img), min_distance=5) + + extractor = BRIEF(descriptor_size=8, sigma=2, mode='uniform') + + extractor.extract(img, keypoints[:8]) + + expected = np.array([[ True, False, True, False, False, True, False, False], + [False, True, False, False, True, True, True, True], + [ True, False, False, False, False, False, False, False], + [False, True, True, False, False, False, True, False], + [False, False, False, False, False, False, True, False], + [False, True, False, False, True, False, False, False], + [False, False, True, True, False, False, True, True], + [ True, True, False, False, False, False, False, False]], dtype=bool) + + assert_array_equal(extractor.descriptors, expected) + + +def test_unsupported_mode(): + assert_raises(ValueError, BRIEF, mode='foobar') + + +def test_border(): + img = np.zeros((100, 100)) + keypoints = np.array([[1, 1], [20, 20], [50, 50], [80, 80]]) + + extractor = BRIEF(patch_size=41) + extractor.extract(img, keypoints) + + assert extractor.descriptors.shape[0] == 3 + assert_array_equal(extractor.mask, (False, True, True, True)) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() diff --git a/skimage/feature/tests/test_censure.py b/skimage/feature/tests/test_censure.py new file mode 100644 index 00000000..f1d88f3d --- /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 CENSURE + + +img = moon() + + +def test_keypoints_censure_color_image_unsupported_error(): + """Censure keypoints can be extracted from gray-scale images only.""" + assert_raises(ValueError, CENSURE().detect, np.zeros((20, 20, 3))) + + +def test_keypoints_censure_mode_validity_error(): + """Mode argument in keypoints_censure can be either DoB, Octagon or + STAR.""" + assert_raises(ValueError, CENSURE, 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.""" + assert_raises(ValueError, CENSURE, 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.""" + detector = CENSURE() + detector.detect(img) + expected_keypoints = np.array([[ 21, 497], + [ 36, 46], + [119, 350], + [185, 177], + [287, 250], + [357, 239], + [463, 116], + [464, 132], + [467, 260]]) + expected_scales = np.array([3, 4, 4, 2, 2, 3, 2, 2, 2]) + + assert_array_equal(expected_keypoints, detector.keypoints) + assert_array_equal(expected_scales, detector.scales) + + +def test_keypoints_censure_moon_image_octagon(): + """Verify the actual Censure keypoints and their corresponding scale with + the expected values for Octagon filter.""" + + detector = CENSURE(mode='octagon') + detector.detect(img) + expected_keypoints = np.array([[ 21, 496], + [ 35, 46], + [287, 250], + [356, 239], + [463, 116]]) + + expected_scales = np.array([3, 4, 2, 2, 2]) + + assert_array_equal(expected_keypoints, detector.keypoints) + assert_array_equal(expected_scales, detector.scales) + + +def test_keypoints_censure_moon_image_star(): + """Verify the actual Censure keypoints and their corresponding scale with + the expected values for STAR filter.""" + detector = CENSURE(mode='star') + detector.detect(img) + expected_keypoints = np.array([[ 21, 497], + [ 36, 46], + [117, 356], + [185, 177], + [260, 227], + [287, 250], + [357, 239], + [451, 281], + [463, 116], + [467, 260]]) + + expected_scales = np.array([3, 3, 6, 2, 3, 2, 3, 5, 2, 2]) + + assert_array_equal(expected_keypoints, detector.keypoints) + assert_array_equal(expected_scales, detector.scales) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index 7ff34796..21070779 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -1,12 +1,94 @@ import numpy as np -from numpy.testing import assert_array_equal, assert_almost_equal +from numpy.testing import (assert_array_equal, assert_raises, + assert_almost_equal) from skimage import data from skimage import img_as_float +from skimage.color import rgb2gray +from skimage.morphology import octagon from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi, corner_subpix, peak_local_max, corner_peaks, - corner_kitchen_rosenfeld, corner_foerstner) + corner_kitchen_rosenfeld, corner_foerstner, + corner_fast, corner_orientations, + structure_tensor, structure_tensor_eigvals, + hessian_matrix, hessian_matrix_eigvals) + + +def test_structure_tensor(): + square = np.zeros((5, 5)) + square[2, 2] = 1 + Axx, Axy, Ayy = structure_tensor(square, sigma=0.1) + assert_array_equal(Axx, np.array([[ 0, 0, 0, 0, 0], + [ 0, 1, 0, 1, 0], + [ 0, 4, 0, 4, 0], + [ 0, 1, 0, 1, 0], + [ 0, 0, 0, 0, 0]])) + assert_array_equal(Axy, np.array([[ 0, 0, 0, 0, 0], + [ 0, 1, 0, -1, 0], + [ 0, 0, 0, -0, 0], + [ 0, -1, -0, 1, 0], + [ 0, 0, 0, 0, 0]])) + assert_array_equal(Ayy, np.array([[ 0, 0, 0, 0, 0], + [ 0, 1, 4, 1, 0], + [ 0, 0, 0, 0, 0], + [ 0, 1, 4, 1, 0], + [ 0, 0, 0, 0, 0]])) + + +def test_hessian_matrix(): + square = np.zeros((5, 5)) + square[2, 2] = 1 + Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) + assert_array_equal(Hxx, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_array_equal(Hxy, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_array_equal(Hyy, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + + +def test_structure_tensor_eigvals(): + square = np.zeros((5, 5)) + square[2, 2] = 1 + Axx, Axy, Ayy = structure_tensor(square, sigma=0.1) + l1, l2 = structure_tensor_eigvals(Axx, Axy, Ayy) + assert_array_equal(l1, np.array([[0, 0, 0, 0, 0], + [0, 2, 4, 2, 0], + [0, 4, 0, 4, 0], + [0, 2, 4, 2, 0], + [0, 0, 0, 0, 0]])) + assert_array_equal(l2, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + + +def test_hessian_matrix_eigvals(): + square = np.zeros((5, 5)) + square[2, 2] = 1 + Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) + l1, l2 = hessian_matrix_eigvals(Hxx, Hxy, Hyy) + assert_array_equal(l1, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_array_equal(l2, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) def test_square_image(): @@ -100,8 +182,8 @@ def test_rotated_lena(): def test_subpix(): img = np.zeros((50, 50)) - img[:25,:25] = 255 - img[25:,25:] = 255 + img[:25, :25] = 255 + img[25:, 25:] = 255 corner = peak_local_max(corner_harris(img), num_peaks=1) subpix = corner_subpix(img, corner) assert_array_equal(subpix[0], (24.5, 24.5)) @@ -128,7 +210,7 @@ def test_num_peaks(): peak_local_max returns exactly the right amount of peaks. Test is run on Lena in order to produce a sufficient number of corners""" - lena_corners = corner_harris(data.lena()) + lena_corners = corner_harris(rgb2gray(data.lena())) for i in range(20): n = np.random.random_integers(20) @@ -166,6 +248,59 @@ def test_blank_image_nans(): assert np.all(np.isfinite(response)) +def test_corner_fast_image_unsupported_error(): + img = np.zeros((20, 20, 3)) + assert_raises(ValueError, corner_fast, img) + + +def test_corner_fast_lena(): + img = rgb2gray(data.lena()) + expected = np.array([[ 67, 157], + [204, 261], + [247, 146], + [269, 111], + [318, 158], + [386, 73], + [413, 70], + [435, 180], + [455, 177], + [461, 160]]) + actual = corner_peaks(corner_fast(img, 12, 0.3)) + assert_array_equal(actual, expected) + + +def test_corner_orientations_image_unsupported_error(): + img = np.zeros((20, 20, 3)) + assert_raises(ValueError, corner_orientations, img, + np.asarray([[7, 7]]), np.ones((3, 3))) + + +def test_corner_orientations_even_shape_error(): + img = np.zeros((20, 20)) + assert_raises(ValueError, corner_orientations, img, + np.asarray([[7, 7]]), np.ones((4, 4))) + + +def test_corner_orientations_lena(): + img = rgb2gray(data.lena()) + corners = corner_peaks(corner_fast(img, 11, 0.35)) + expected = np.array([-1.9195897 , -3.03159624, -1.05991162, -2.89573739, + -2.61607644, 2.98660159]) + actual = corner_orientations(img, corners, octagon(3, 2)) + assert_almost_equal(actual, expected) + + +def test_corner_orientations_square(): + square = np.zeros((12, 12)) + square[3:9, 3:9] = 1 + corners = corner_peaks(corner_fast(square, 9), min_distance=1) + actual_orientations = corner_orientations(square, corners, octagon(3, 2)) + actual_orientations_degrees = np.rad2deg(actual_orientations) + expected_orientations_degree = np.array([ 45., 135., -45., -135.]) + assert_array_equal(actual_orientations_degrees, + expected_orientations_degree) + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() diff --git a/skimage/feature/tests/test_match.py b/skimage/feature/tests/test_match.py new file mode 100644 index 00000000..f9496c0c --- /dev/null +++ b/skimage/feature/tests/test_match.py @@ -0,0 +1,96 @@ +import numpy as np +from numpy.testing import assert_equal, assert_raises +from skimage import data +from skimage import transform as tf +from skimage.color import rgb2gray +from skimage.feature import (BRIEF, match_descriptors, + corner_peaks, corner_harris) + + +def test_binary_descriptors_unequal_descriptor_sizes_error(): + """Sizes of descriptors of keypoints to be matched should be equal.""" + descs1 = np.array([[True, True, False, True], + [False, True, False, True]]) + descs2 = np.array([[True, False, False, True, False], + [False, True, True, True, False]]) + assert_raises(ValueError, match_descriptors, descs1, descs2) + + +def test_binary_descriptors(): + descs1 = np.array([[True, True, False, True, True], + [False, True, False, True, True]]) + descs2 = np.array([[True, False, False, True, False], + [False, False, True, True, True]]) + matches = match_descriptors(descs1, descs2) + assert_equal(matches, [[0, 0], [1, 1]]) + + +def test_binary_descriptors_lena_rotation_crosscheck_false(): + """Verify matched keypoints and their corresponding masks results between + lena image and its rotated version with the expected keypoint pairs with + cross_check disabled.""" + img = data.lena() + img = rgb2gray(img) + tform = tf.SimilarityTransform(scale=1, rotation=0.15, translation=(0, 0)) + rotated_img = tf.warp(img, tform) + + extractor = BRIEF(descriptor_size=512) + + keypoints1 = corner_peaks(corner_harris(img), min_distance=5) + extractor.extract(img, keypoints1) + descriptors1 = extractor.descriptors + + keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5) + extractor.extract(rotated_img, keypoints2) + descriptors2 = extractor.descriptors + + matches = match_descriptors(descriptors1, descriptors2, threshold=0.13, + cross_check=False) + + exp_matches1 = np.array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, + 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, + 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, + 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]) + exp_matches2 = np.array([33, 0, 35, 7, 1, 35, 3, 2, 3, 6, 4, 9, + 11, 10, 28, 7, 8, 5, 31, 14, 13, 15, 21, 16, + 16, 13, 17, 18, 19, 21, 22, 23, 0, 24, 1, 24, + 23, 0, 26, 27, 25, 34, 28, 14, 29, 30, 21]) + assert_equal(matches[:, 0], exp_matches1) + assert_equal(matches[:, 1], exp_matches2) + + +def test_binary_descriptors_lena_rotation_crosscheck_true(): + """Verify matched keypoints and their corresponding masks results between + lena image and its rotated version with the expected keypoint pairs with + cross_check enabled.""" + img = data.lena() + img = rgb2gray(img) + tform = tf.SimilarityTransform(scale=1, rotation=0.15, translation=(0, 0)) + rotated_img = tf.warp(img, tform) + + extractor = BRIEF(descriptor_size=512) + + keypoints1 = corner_peaks(corner_harris(img), min_distance=5) + extractor.extract(img, keypoints1) + descriptors1 = extractor.descriptors + + keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5) + extractor.extract(rotated_img, keypoints2) + descriptors2 = extractor.descriptors + + matches = match_descriptors(descriptors1, descriptors2, threshold=0.13, + cross_check=True) + + exp_matches1 = np.array([ 0, 1, 2, 4, 6, 7, 9, 10, 11, 12, 13, 15, + 16, 17, 19, 20, 21, 24, 26, 27, 28, 29, 30, 35, + 36, 38, 39, 40, 42, 44, 45]) + exp_matches2 = np.array([33, 0, 35, 1, 3, 2, 6, 4, 9, 11, 10, 7, + 8, 5, 14, 13, 15, 16, 17, 18, 19, 21, 22, 24, + 23, 26, 27, 25, 28, 29, 30]) + assert_equal(matches[:, 0], exp_matches1) + assert_equal(matches[:, 1], exp_matches2) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() diff --git a/skimage/feature/tests/test_orb.py b/skimage/feature/tests/test_orb.py new file mode 100644 index 00000000..30394d07 --- /dev/null +++ b/skimage/feature/tests/test_orb.py @@ -0,0 +1,115 @@ +import numpy as np +from numpy.testing import assert_array_equal, assert_almost_equal +from skimage.feature import ORB +from skimage.data import lena +from skimage.color import rgb2gray + + +img = rgb2gray(lena()) + + +def test_keypoints_orb_desired_no_of_keypoints(): + detector_extractor = ORB(n_keypoints=10, fast_n=12, fast_threshold=0.20) + detector_extractor.detect(img) + + exp_rows = np.array([ 435. , 435.6 , 376. , 455. , 434.88, 269. , + 375.6 , 310.8 , 413. , 311.04]) + exp_cols = np.array([ 180. , 180. , 156. , 176. , 180. , 111. , + 156. , 172.8, 70. , 172.8]) + + exp_scales = np.array([ 1. , 1.2 , 1. , 1. , 1.44 , 1. , + 1.2 , 1.2 , 1. , 1.728]) + + exp_orientations = np.array([-175.64733392, -167.94842949, -148.98350192, + -142.03599837, -176.08535837, -53.08162354, + -150.89208271, 97.7693776 , -173.4479964 , + 38.66312042]) + exp_response = np.array([ 0.96770745, 0.81027306, 0.72376257, + 0.5626413 , 0.5097993 , 0.44351774, + 0.39154173, 0.39084861, 0.39063076, + 0.37602487]) + + assert_almost_equal(exp_rows, detector_extractor.keypoints[:, 0]) + assert_almost_equal(exp_cols, detector_extractor.keypoints[:, 1]) + assert_almost_equal(exp_scales, detector_extractor.scales) + assert_almost_equal(exp_response, detector_extractor.responses) + assert_almost_equal(exp_orientations, + np.rad2deg(detector_extractor.orientations), 5) + + detector_extractor.detect_and_extract(img) + assert_almost_equal(exp_rows, detector_extractor.keypoints[:, 0]) + assert_almost_equal(exp_cols, detector_extractor.keypoints[:, 1]) + + +def test_keypoints_orb_less_than_desired_no_of_keypoints(): + img = rgb2gray(lena()) + detector_extractor = ORB(n_keypoints=15, fast_n=12, + fast_threshold=0.33, downscale=2, n_scales=2) + detector_extractor.detect(img) + + exp_rows = np.array([ 67., 247., 269., 413., 435., 230., 264., + 330., 372.]) + exp_cols = np.array([ 157., 146., 111., 70., 180., 136., 336., + 148., 156.]) + + exp_scales = np.array([ 1., 1., 1., 1., 1., 2., 2., 2., 2.]) + + exp_orientations = np.array([-105.76503839, -96.28973044, -53.08162354, + -173.4479964 , -175.64733392, -106.07927215, + -163.40016243, 75.80865813, -154.73195911]) + + exp_response = np.array([ 0.13197835, 0.24931321, 0.44351774, + 0.39063076, 0.96770745, 0.04935129, + 0.21431068, 0.15826555, 0.42403573]) + + assert_almost_equal(exp_rows, detector_extractor.keypoints[:, 0]) + assert_almost_equal(exp_cols, detector_extractor.keypoints[:, 1]) + assert_almost_equal(exp_scales, detector_extractor.scales) + assert_almost_equal(exp_response, detector_extractor.responses) + assert_almost_equal(exp_orientations, + np.rad2deg(detector_extractor.orientations), 5) + + detector_extractor.detect_and_extract(img) + assert_almost_equal(exp_rows, detector_extractor.keypoints[:, 0]) + assert_almost_equal(exp_cols, detector_extractor.keypoints[:, 1]) + + +def test_descriptor_orb(): + detector_extractor = ORB(fast_n=12, fast_threshold=0.20) + + exp_descriptors = np.array([[ True, False, True, True, False, False, False, False, False, False], + [False, False, True, True, False, True, True, False, True, True], + [ True, False, False, False, True, False, True, True, True, False], + [ True, False, False, True, False, True, True, False, False, False], + [False, True, True, True, False, False, False, True, True, False], + [False, False, False, False, False, True, False, True, True, True], + [False, True, True, True, True, False, False, True, False, True], + [ True, True, True, False, True, True, True, True, False, False], + [ True, True, False, True, True, True, True, False, False, False], + [ True, False, False, False, False, True, False, False, True, True], + [ True, False, False, False, True, True, True, False, False, False], + [False, False, True, False, True, False, False, True, False, False], + [False, False, True, True, False, False, False, False, False, True], + [ True, True, False, False, False, True, True, True, True, True], + [ True, True, True, False, False, True, False, True, True, False], + [False, True, True, False, False, True, True, True, True, True], + [ True, True, True, False, False, False, False, True, True, True], + [False, False, False, False, True, False, False, True, True, False], + [False, True, False, False, True, False, False, False, True, True], + [ True, False, True, False, False, False, True, True, False, False]], dtype=bool) + + detector_extractor.detect(img) + detector_extractor.extract(img, detector_extractor.keypoints, + detector_extractor.scales, + detector_extractor.orientations) + assert_array_equal(exp_descriptors, + detector_extractor.descriptors[100:120, 10:20]) + + detector_extractor.detect_and_extract(img) + assert_array_equal(exp_descriptors, + detector_extractor.descriptors[100:120, 10:20]) + + +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 index 6e25f51a..2b601b94 100644 --- a/skimage/feature/tests/test_util.py +++ b/skimage/feature/tests/test_util.py @@ -1,30 +1,71 @@ import numpy as np -from numpy.testing import assert_array_equal -from skimage.feature.util import pairwise_hamming_distance +import matplotlib.pyplot as plt +from numpy.testing import assert_equal, assert_raises + +from skimage.feature.util import (FeatureDetector, DescriptorExtractor, + _prepare_grayscale_input_2D, + _mask_border_keypoints, plot_matches) -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_feature_detector(): + assert_raises(NotImplementedError, FeatureDetector().detect, None) -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) +def test_descriptor_extractor(): + assert_raises(NotImplementedError, DescriptorExtractor().extract, + None, None) + + +def test_prepare_grayscale_input_2D(): + assert_raises(ValueError, _prepare_grayscale_input_2D, np.zeros((3, 3, 3))) + assert_raises(ValueError, _prepare_grayscale_input_2D, np.zeros((3, 1))) + assert_raises(ValueError, _prepare_grayscale_input_2D, np.zeros((3, 1, 1))) + img = _prepare_grayscale_input_2D(np.zeros((3, 3))) + img = _prepare_grayscale_input_2D(np.zeros((3, 3, 1))) + img = _prepare_grayscale_input_2D(np.zeros((1, 3, 3))) + + +def test_mask_border_keypoints(): + keypoints = np.array([[0, 0], [1, 1], [2, 2], [3, 3], [4, 4]]) + assert_equal(_mask_border_keypoints((10, 10), keypoints, 0), + [1, 1, 1, 1, 1]) + assert_equal(_mask_border_keypoints((10, 10), keypoints, 2), + [0, 0, 1, 1, 1]) + assert_equal(_mask_border_keypoints((4, 4), keypoints, 2), + [0, 0, 1, 0, 0]) + assert_equal(_mask_border_keypoints((10, 10), keypoints, 5), + [0, 0, 0, 0, 0]) + assert_equal(_mask_border_keypoints((10, 10), keypoints, 4), + [0, 0, 0, 0, 1]) + + +def test_plot_matches(): + fig, ax = plt.subplots(nrows=1, ncols=1) + + shapes = (((10, 10), (10, 10)), + ((10, 10), (12, 10)), + ((10, 10), (10, 12)), + ((10, 10), (12, 12)), + ((12, 10), (10, 10)), + ((10, 12), (10, 10)), + ((12, 12), (10, 10))) + + keypoints1 = 10 * np.random.rand(10, 2) + keypoints2 = 10 * np.random.rand(10, 2) + idxs1 = np.random.randint(10, size=10) + idxs2 = np.random.randint(10, size=10) + matches = np.column_stack((idxs1, idxs2)) + + for shape1, shape2 in shapes: + img1 = np.zeros(shape1) + img2 = np.zeros(shape2) + plot_matches(ax, img1, img2, keypoints1, keypoints2, matches) + plot_matches(ax, img1, img2, keypoints1, keypoints2, matches, + only_matches=True) + plot_matches(ax, img1, img2, keypoints1, keypoints2, matches, + keypoints_color='r') + plot_matches(ax, img1, img2, keypoints1, keypoints2, matches, + matches_color='r') if __name__ == '__main__': diff --git a/skimage/feature/util.py b/skimage/feature/util.py index a5267d44..8ee2baf8 100644 --- a/skimage/feature/util.py +++ b/skimage/feature/util.py @@ -1,38 +1,161 @@ +import numpy as np + +from skimage.util import img_as_float -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] +class FeatureDetector(object): - keypoints_filtering_mask = ((dist - 1 < keypoints[:, 0]) & - (keypoints[:, 0] < width - dist + 1) & - (dist - 1 < keypoints[:, 1]) & - (keypoints[:, 1] < height - dist + 1)) + def __init__(self): + self.keypoints_ = np.array([]) - return keypoints_filtering_mask + def detect(self, image): + """Detect keypoints in image. + + Parameters + ---------- + image : 2D array + Input image. + + """ + raise NotImplementedError() -def pairwise_hamming_distance(array1, array2): - """**Experimental function**. +class DescriptorExtractor(object): - Calculate hamming dissimilarity measure between two sets of - vectors. + def __init__(self): + self.descriptors_ = np.array([]) + + def extract(self, image, keypoints): + """Extract feature descriptors in image for given keypoints. + + Parameters + ---------- + image : 2D array + Input image. + keypoints : (N, 2) array + Keypoint locations as ``(row, col)``. + + """ + raise NotImplementedError() + + +def plot_matches(ax, image1, image2, keypoints1, keypoints2, matches, + keypoints_color='k', matches_color=None, only_matches=False): + """Plot matched features. Parameters ---------- - array1 : (P1, D) array - P1 vectors of size D. - array2 : (P2, D) array - P2 vectors of size D. + ax : matplotlib.axes.Axes + Matches and image are drawn in this ax. + image1 : (N, M [, 3]) array + First grayscale or color image. + image2 : (N, M [, 3]) array + Second grayscale or color image. + keypoints1 : (K1, 2) array + First keypoint coordinates as ``(row, col)``. + keypoints2 : (K2, 2) array + Second keypoint coordinates as ``(row, col)``. + matches : (Q, 2) array + Indices of corresponding matches in first and second set of + descriptors, where ``matches[:, 0]`` denote the indices in the first + and ``matches[:, 1]`` the indices in the second set of descriptors. + keypoints_color : matplotlib color, optional + Color for keypoint locations. + matches_color : matplotlib color, optional + Color for lines which connect keypoint matches. By default the + color is chosen randomly. + only_matches : bool, optional + Whether to only plot matches and not plot the keypoint locations. + + """ + + image1 = img_as_float(image1) + image2 = img_as_float(image2) + + new_shape1 = list(image1.shape) + new_shape2 = list(image2.shape) + + if image1.shape[0] < image2.shape[0]: + new_shape1[0] = image2.shape[0] + elif image1.shape[0] > image2.shape[0]: + new_shape2[0] = image1.shape[0] + + if image1.shape[1] < image2.shape[1]: + new_shape1[1] = image2.shape[1] + elif image1.shape[1] > image2.shape[1]: + new_shape2[1] = image1.shape[1] + + if new_shape1 != image1.shape: + new_image1 = np.zeros(new_shape1, dtype=image1.dtype) + new_image1[:image1.shape[0], :image1.shape[1]] = image1 + image1 = new_image1 + + if new_shape2 != image2.shape: + new_image2 = np.zeros(new_shape2, dtype=image2.dtype) + new_image2[:image2.shape[0], :image2.shape[1]] = image2 + image2 = new_image2 + + image = np.concatenate([image1, image2], axis=1) + + offset = image1.shape + + if not only_matches: + ax.scatter(keypoints1[:, 1], keypoints1[:, 0], + facecolors='none', edgecolors=keypoints_color) + ax.scatter(keypoints2[:, 1] + offset[1], keypoints2[:, 0], + facecolors='none', edgecolors=keypoints_color) + + ax.imshow(image) + ax.axis((0, 2 * offset[1], offset[0], 0)) + + for i in range(matches.shape[0]): + idx1 = matches[i, 0] + idx2 = matches[i, 1] + + if matches_color is None: + color = np.random.rand(3, 1) + else: + color = matches_color + + ax.plot((keypoints1[idx1, 1], keypoints2[idx2, 1] + offset[1]), + (keypoints1[idx1, 0], keypoints2[idx2, 0]), + '-', color=color) + + +def _prepare_grayscale_input_2D(image): + image = np.squeeze(image) + if image.ndim != 2: + raise ValueError("Only 2-D gray-scale images supported.") + + return img_as_float(image) + + +def _mask_border_keypoints(image_shape, keypoints, distance): + """Mask coordinates that are within certain distance from the image border. + + Parameters + ---------- + image_shape : (2, ) array_like + Shape of the image as ``(rows, cols)``. + keypoints : (N, 2) array + Keypoint coordinates as ``(rows, cols)``. + distance : int + Image border distance. 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. + mask : (N, ) bool array + Mask indicating if pixels are within the image (``True``) or in the + border region of the image (``False``). """ - distance = (array1[:, None] != array2[None]).mean(axis=2) - return distance + + rows = image_shape[0] + cols = image_shape[1] + + mask = (((distance - 1) < keypoints[:, 0]) + & (keypoints[:, 0] < (rows - distance + 1)) + & ((distance - 1) < keypoints[:, 1]) + & (keypoints[:, 1] < (cols - distance + 1))) + + return mask