wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-25 13:30:51 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #668 from ankit-maverick/censure

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