From 0b37611df60bf1da7542324f77e9d52093aab824 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 11 Aug 2013 11:02:24 +0200 Subject: [PATCH] Fix several bugs in DoB method and improve overall code quality --- skimage/feature/censure.py | 76 +++++++++++++++------------ skimage/feature/censure_cy.pyx | 21 ++++---- skimage/feature/tests/test_censure.py | 15 +++--- 3 files changed, 59 insertions(+), 53 deletions(-) diff --git a/skimage/feature/censure.py b/skimage/feature/censure.py index e3264436..583e3d0a 100644 --- a/skimage/feature/censure.py +++ b/skimage/feature/censure.py @@ -24,7 +24,16 @@ def _get_filtered_image(image, n_scales, mode): scales = np.zeros((image.shape[0], image.shape[1], n_scales), dtype=np.double) - if mode == 'DoB': + if mode == 'dob': + + # make scales[:, :, i] contiguous memory block + item_size = scales.itemsize + scales.strides = (item_size * scales.shape[0], + item_size, + item_size * scales.shape[0] * scales.shape[1]) + + integral_img = integral_image(image) + for i in range(n_scales): n = i + 1 @@ -34,33 +43,29 @@ def _get_filtered_image(image, n_scales, mode): inner_weight = (1.0 / (2 * n + 1)**2) outer_weight = (1.0 / (12 * n**2 + 4 * n)) - integral_img = integral_image(image) - - filtered_image = np.zeros(image.shape) - _censure_dob_loop(image, n, integral_img, filtered_image, + _censure_dob_loop(n, integral_img, scales[:, :, i], inner_weight, outer_weight) - scales[:, :, i] = filtered_image - # NOTE : For the Octagon shaped filter, we implemented and evaluated the # slanted integral image based image filtering but the performance was # more or less equal to image filtering using # scipy.ndimage.filters.convolve(). Hence we have decided to use the # later for a much cleaner implementation. - elif mode == 'Octagon': + elif mode == 'octagon': # TODO : Decide the shapes of Octagon filters for scales > 7 for i in range(n_scales): + mo, no = OCTAGON_OUTER_SHAPE[i] + mi, ni = OCTAGON_INNER_SHAPE[i] scales[:, :, i] = convolve(image, - _octagon_filter_kernel(OCTAGON_OUTER_SHAPE[i][0], - OCTAGON_OUTER_SHAPE[i][1], OCTAGON_INNER_SHAPE[i][0], - OCTAGON_INNER_SHAPE[i][1])) - else: + _octagon_filter_kernel(mo, no, mi, ni)) + elif mode == 'star': for i in range(n_scales): + m = STAR_SHAPE[STAR_FILTER_SHAPE[i][0]] + n = STAR_SHAPE[STAR_FILTER_SHAPE[i][1]] scales[:, :, i] = convolve(image, - _star_filter_kernel(STAR_SHAPE[STAR_FILTER_SHAPE[i][0]], - STAR_SHAPE[STAR_FILTER_SHAPE[i][1]])) + _star_filter_kernel(m, n)) return scales @@ -115,7 +120,7 @@ def _star(a): def _star_filter_kernel(m, n): c = m + m // 2 - n - n // 2 outer_star = _star(m) - inner_star = np.zeros((outer_star.shape)) + inner_star = np.zeros_like(outer_star) inner_star[c: -c, c: -c] = _star(n) outer_weight = 1.0 / (np.sum(outer_star - inner_star)) inner_weight = 1.0 / np.sum(inner_star) @@ -128,7 +133,6 @@ 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 - return feature_mask def censure_keypoints(image, n_scales=7, mode='DoB', non_max_threshold=0.15, @@ -141,19 +145,15 @@ def censure_keypoints(image, n_scales=7, mode='DoB', non_max_threshold=0.15, ---------- image : 2D ndarray Input image. - n_scales : positive integer Number of scales to extract keypoints from. The keypoints will be extracted from all the scales except the first and the last. - mode : ('DoB', 'Octagon', 'STAR') Type of bilevel filter used to get the scales of input image. Possible values are 'DoB', 'Octagon' and 'STAR'. - non_max_threshold : float Threshold value used to suppress maximas and minimas with a weak magnitude response obtained after Non-Maximal Suppression. - line_threshold : float Threshold for rejecting interest points which have ratio of principal curvatures greater than this value. @@ -162,8 +162,7 @@ def censure_keypoints(image, n_scales=7, mode='DoB', non_max_threshold=0.15, ------- keypoints : (N, 2) array Location of the extracted keypoints in the (row, col) format. - - scale : (N, 1) array + scales : (N, 1) array The corresponding scale of the N extracted keypoints. References @@ -183,10 +182,15 @@ def censure_keypoints(image, n_scales=7, mode='DoB', non_max_threshold=0.15, image = np.squeeze(image) if image.ndim != 2: raise ValueError("Only 2-D gray-scale images supported.") - image = img_as_float(image) + image = img_as_float(image) image = np.ascontiguousarray(image) + mode = mode.lower() + + if mode not in ('dob', 'octagon', 'star'): + raise ValueError('Mode must be one of "DoB", "Octagon", "STAR".') + # Generating all the scales filter_response = _get_filtered_image(image, n_scales, mode) @@ -199,29 +203,33 @@ def censure_keypoints(image, n_scales=7, mode='DoB', non_max_threshold=0.15, feature_mask[filter_response < non_max_threshold] = False for i in range(1, n_scales - 1): - # sigma = (window_size - 1) / 6.0 + # sigma = (window_size - 1) / 6.0, so the window covers > 99% of the + # kernel's distribution # window_size = 7 + 2 * i # Hence sigma = 1 + i / 3.0 - feature_mask[:, :, i] = _suppress_lines(feature_mask[:, :, i], image, - (1 + i / 3.0), line_threshold) + _suppress_lines(feature_mask[:, :, i], image, + (1 + i / 3.0), line_threshold) rows, cols, scales = np.nonzero(feature_mask[..., 1:n_scales - 1]) keypoints = np.column_stack([rows, cols]) scales = scales + 2 - if mode == 'DoB': + if mode == 'dob': return keypoints, scales cumulative_mask = np.zeros(keypoints.shape[0], dtype=np.bool) - if mode == 'Octagon': + if mode == 'octagon': for i in range(2, n_scales): - c = (OCTAGON_OUTER_SHAPE[i - 1][0] - 1) // 2 + OCTAGON_OUTER_SHAPE[i - 1][1] - cumulative_mask = cumulative_mask | (_mask_border_keypoints(image, keypoints, c) & (scales == i)) - - elif mode == 'STAR': + 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(2, n_scales): - c = STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] + STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] // 2 - cumulative_mask = cumulative_mask | (_mask_border_keypoints(image, keypoints, c) & (scales == i)) + c = STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] \ + + STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] // 2 + cumulative_mask |= _mask_border_keypoints(image, keypoints, c) \ + & (scales == i) return keypoints[cumulative_mask], scales[cumulative_mask] diff --git a/skimage/feature/censure_cy.pyx b/skimage/feature/censure_cy.pyx index 93c7c142..cfd1260f 100644 --- a/skimage/feature/censure_cy.pyx +++ b/skimage/feature/censure_cy.pyx @@ -3,29 +3,28 @@ #cython: nonecheck=False #cython: wraparound=False -cimport numpy as cnp -import numpy as np - -def _censure_dob_loop(double[:, ::1] image, Py_ssize_t n, +def _censure_dob_loop(Py_ssize_t n, double[:, ::1] integral_img, double[:, ::1] filtered_image, double inner_weight, double outer_weight): cdef Py_ssize_t i, j cdef double inner, outer + cdef Py_ssize_t n2 = 2 * n + cdef double total_weight = inner_weight + outer_weight - for i in range(2 * n, image.shape[0] - 2 * n): - for j in range(2 * n, image.shape[1] - 2 * n): + for i in range(n2, integral_img.shape[0] - n2): + for j in range(n2, integral_img.shape[1] - n2): inner = (integral_img[i + n, j + n] + integral_img[i - n - 1, j - n - 1] - integral_img[i + n, j - n - 1] - integral_img[i - n - 1, j + n]) - outer = (integral_img[i + 2 * n, j + 2 * n] - + integral_img[i - 2 * n - 1, j - 2 * n - 1] - - integral_img[i + 2 * n, j - 2 * n - 1] - - integral_img[i - 2 * n - 1, j + 2 * n]) + outer = (integral_img[i + n2, j + n2] + + integral_img[i - n2 - 1, j - n2 - 1] + - integral_img[i + n2, j - n2 - 1] + - integral_img[i - n2 - 1, j + n2]) filtered_image[i, j] = (outer_weight * outer - - (inner_weight + outer_weight) * inner) + - total_weight * inner) diff --git a/skimage/feature/tests/test_censure.py b/skimage/feature/tests/test_censure.py index 41d1886e..14e97a39 100644 --- a/skimage/feature/tests/test_censure.py +++ b/skimage/feature/tests/test_censure.py @@ -1,6 +1,7 @@ import numpy as np from numpy.testing import assert_array_equal, assert_raises from skimage.data import moon +from skimage.util import img_as_ubyte from skimage.feature import censure_keypoints @@ -15,10 +16,7 @@ def test_censure_keypoints_moon_image_DoB(): the expected values for DoB filter.""" img = moon() actual_kp_DoB, actual_scale = censure_keypoints(img, 7, 'DoB', 0.15) - expected_kp_DoB = np.array([[ 4, 507], - [ 8, 503], - [ 12, 499], - [ 21, 497], + expected_kp_DoB = np.array([[ 21, 497], [ 36, 46], [119, 350], [185, 177], @@ -27,7 +25,7 @@ def test_censure_keypoints_moon_image_DoB(): [463, 116], [464, 132], [467, 260]]) - expected_scale = np.array([2, 4, 6, 3, 4, 4, 2, 2, 3, 2, 2, 2]) + 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) @@ -37,14 +35,15 @@ def test_censure_keypoints_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 = censure_keypoints(img, 7, 'Octagon', 0.15) + actual_kp_Octagon, actual_scale = censure_keypoints(img, 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], dtype=np.intp) + 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) @@ -66,7 +65,7 @@ def test_censure_keypoints_moon_image_STAR(): [463, 116], [467, 260]]) - expected_scale = np.array([3, 3, 6, 2, 3, 2, 3, 5, 2, 2], dtype=np.intp) + 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)