diff --git a/bento.info b/bento.info index 45f57a36..f1cbabc1 100644 --- a/bento.info +++ b/bento.info @@ -151,6 +151,9 @@ Library: Extension: skimage.restoration._denoise_cy Sources: skimage/restoration/_denoise_cy.pyx + Extension: skimage.feature._hessian_det_appx + Sources: + skimage/exposure/_hessian_det_appx.pyx Executable: skivi Module: skimage.scripts.skivi diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 54b7bac3..c46fde01 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -6,7 +6,7 @@ from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi, corner_foerstner, corner_subpix, corner_peaks, corner_fast, structure_tensor, structure_tensor_eigvals, hessian_matrix, - hessian_matrix_eigvals) + hessian_matrix_eigvals, hessian_matrix_det) from .corner_cy import corner_moravec, corner_orientations from .template import match_template from .brief import BRIEF @@ -14,7 +14,7 @@ from .censure import CENSURE from .orb import ORB from .match import match_descriptors from .util import plot_matches -from .blob import blob_dog, blob_log +from .blob import blob_dog, blob_log, blob_doh __all__ = ['daisy', @@ -26,6 +26,7 @@ __all__ = ['daisy', 'structure_tensor', 'structure_tensor_eigvals', 'hessian_matrix', + 'hessian_matrix_det', 'hessian_matrix_eigvals', 'corner_kitchen_rosenfeld', 'corner_harris', @@ -43,4 +44,5 @@ __all__ = ['daisy', 'match_descriptors', 'plot_matches', 'blob_dog', + 'blob_doh', 'blob_log'] diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx new file mode 100644 index 00000000..0c54cba7 --- /dev/null +++ b/skimage/feature/_hessian_det_appx.pyx @@ -0,0 +1,153 @@ +# cython: cdivision=True +# cython: boundscheck=False +# cython: nonecheck=False +# cython: wraparound=False +import numpy as np +cimport numpy as cnp + + +cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high): + """Clips coordinate between high and low. + + This method was created so that `hessian_det_appx` does not have to make + a Python call. + + Parameters + ---------- + x : int + Coordinate to be clipped. + low : int + The lower bound. + high : int + The higher bound. + + Returns + ------- + x : int + `x` clipped between `high` and `low`. + """ + + if(x > high): + return high + if(x < low): + return low + return x + + +cdef inline cnp.double_t _integ( + cnp.double_t[:, ::1] img, Py_ssize_t r, Py_ssize_t c, + Py_ssize_t rl, Py_ssize_t cl): + """Integrate over the integral image in the given window + + This method was created so that `hessian_det_appx` does not have to make + a Python call. + + Parameters + ---------- + img : array + The integral image over which to integrate. + r : int + The row number of the top left corner. + c : int + The column number of the top left corner. + rl : int + The number of rows over which to integrate. + cl : int + The number of columns over which to integrate. + + Returns + ------- + ans : int + The integral over the given window. + """ + + r = _clip(r, 0, img.shape[0] - 1) + c = _clip(c, 0, img.shape[1] - 1) + + r2 = _clip(r + rl, 0, img.shape[0] - 1) + c2 = _clip(c + cl, 0, img.shape[1] - 1) + + cdef cnp.double_t ans = img[r, c] + img[r2, c2] - img[r, c2] - img[r2, c] + + if (ans < 0): + return 0 + return ans + + +def _hessian_matrix_det(cnp.double_t[:, ::1] img, double sigma): + """Computes the approximate Hessian Determinant over an image. + + This method uses box filters over integral images to compute the + approximate Hessian Determinant as described in [1]_. + + Parameters + ---------- + img : array + The integral image over which to compute Hessian Determinant. + sigma : float + Standard deviation used for the Gaussian kernel, used for the Hessian + matrix + + Returns + ------- + out : array + The array of the Determinant of Hessians. + + References + ---------- + .. [1] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool, + "SURF: Speeded Up Robust Features" + ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf + + Notes + ----- + The running time of this method only depends on size of the image. It is + independent of `sigma` as one would expect. The downside is that the + result for `sigma` less than `3` is not accurate, i.e., not similar to + the result obtained if someone computed the Hessian and took it's + determinant. + """ + + cdef Py_ssize_t size = int(3 * sigma) + cdef Py_ssize_t height = img.shape[0] + cdef Py_ssize_t width = img.shape[1] + cdef Py_ssize_t r, c + cdef Py_ssize_t s2 = (size - 1) / 2 + cdef Py_ssize_t s3 = size / 3 + cdef Py_ssize_t l = size / 3 + cdef Py_ssize_t w = size + cdef Py_ssize_t b = (size - 1) / 2 + cdef cnp.double_t mid, side, tl, tr, bl, br + cdef cnp.double_t[:, ::1] out = np.zeros_like(img, dtype=np.double) + cdef cnp.double_t w_i = 1.0 / size / size + + cdef float dxx, dyy, dxy + + if not size % 2: + size += 1 + + for r in range(height): + for c in range(width): + tl = _integ(img, r - s3, c - s3, s3, s3) # top left + br = _integ(img, r + 1, c + 1, s3, s3) # bottom right + bl = _integ(img, r - s3, c + 1, s3, s3) # bottom left + tr = _integ(img, r + 1, c - s3, s3, s3) # top right + + dxy = bl + tr - tl - br + dxy = -dxy * w_i + + mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w) # middle box + side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) # sides + + dxx = mid - 3 * side + dxx = -dxx * w_i + + mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1) + side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) + + dyy = mid - 3 * side + dyy = -dyy * w_i + + out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy)) + + return out diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py index 1b465c5f..134e7026 100644 --- a/skimage/feature/blob.py +++ b/skimage/feature/blob.py @@ -1,3 +1,4 @@ + import numpy as np from scipy.ndimage.filters import gaussian_filter, gaussian_laplace import itertools as itt @@ -6,6 +7,8 @@ from math import sqrt, hypot, log from numpy import arccos from skimage.util import img_as_float from .peak import peak_local_max +from ._hessian_det_appx import _hessian_matrix_det +from skimage.transform import integral_image # This basic blob detection algorithm is based on: @@ -33,7 +36,6 @@ def _blob_overlap(blob1, blob2): ------- f : float Fraction of overlapped area. - """ root2 = sqrt(2) @@ -78,7 +80,6 @@ def _prune_blobs(blobs_array, overlap): ------- A : ndarray `array` with overlapping blobs removed. - """ # iterating again might eliminate more blobs, but one iteration suffices @@ -285,8 +286,8 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2, else: sigma_list = np.linspace(min_sigma, max_sigma, num_sigma) - #computing gaussian laplace - #s**2 provides scale invariance + # computing gaussian laplace + # s**2 provides scale invariance gl_images = [-gaussian_laplace(image, s) * s ** 2 for s in sigma_list] image_cube = np.dstack(gl_images) @@ -298,3 +299,112 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2, # Convert the last index to its corresponding scale value local_maxima[:, 2] = sigma_list[local_maxima[:, 2]] return _prune_blobs(local_maxima, overlap) + + +def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01, + overlap=.5, log_scale=False): + """Finds blobs in the given grayscale image. + + Blobs are found using the Determinant of Hessian method [1]_. For each blob + found, the method returns its coordinates and the standard deviation + of the Gaussian Kernel used for the Hessian matrix whose determinant + detected the blob. Determinant of Hessians is approximated using [2]_. + + Parameters + ---------- + image : ndarray + Input grayscale image.Blobs can either be light on dark or vice versa. + min_sigma : float, optional + The minimum standard deviation for Gaussian Kernel used to compute + Hessian matrix. Keep this low to detect smaller blobs. + max_sigma : float, optional + The maximum standard deviation for Gaussian Kernel used to compute + Hessian matrix. Keep this high to detect larger blobs. + num_sigma : int, optional + The number of intermediate values of standard deviations to consider + between `min_sigma` and `max_sigma`. + threshold : float, optional. + The absolute lower bound for scale space maxima. Local maxima smaller + than thresh are ignored. Reduce this to detect less prominent blobs. + overlap : float, optional + A value between 0 and 1. If the area of two blobs overlaps by a + fraction greater than `threshold`, the smaller blob is eliminated. + log_scale : bool, optional + If set intermediate values of standard deviations are interpolated + using a logarithmic scale to the base `10`. If not, linear + interpolation is used. + + Returns + ------- + A : (n, 3) ndarray + A 2d array with each row representing 3 values, ``(y,x,sigma)`` + where ``(y,x)`` are coordinates of the blob and ``sigma`` is the + standard deviation of the Gaussian kernel of the Hessian Matrix whose + determinant detected the blob. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian + + .. [2] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool, + "SURF: Speeded Up Robust Features" + ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf + + Examples + -------- + >>> from skimage import data, feature + >>> img = data.coins() + >>> feature.blob_doh(img) + array([[121, 271, 30], + [123, 44, 23], + [123, 205, 20], + [124, 336, 20], + [126, 101, 20], + [126, 153, 20], + [156, 302, 30], + [185, 348, 30], + [192, 212, 23], + [193, 275, 23], + [195, 100, 23], + [197, 44, 20], + [197, 153, 20], + [260, 173, 30], + [262, 243, 23], + [265, 113, 23], + [270, 363, 30]]) + + + Notes + ----- + The radius of each blob is approximately `sigma`. + Computation of Determinant of Hessians is independent of the standard + deviation. Therefore detecting larger blobs won't take more time. In + methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation + of Gaussians for larger `sigma` takes more time. The downside is that + this method can't be used for detecting blobs of radius less than `3px` + due to the box filters used in the approximation of Hessian Determinant. + """ + + if image.ndim != 2: + raise ValueError("'image' must be grayscale ") + + image = img_as_float(image) + image = integral_image(image) + + if log_scale: + start, stop = log(min_sigma, 10), log(max_sigma, 10) + sigma_list = np.logspace(start, stop, num_sigma) + else: + sigma_list = np.linspace(min_sigma, max_sigma, num_sigma) + + hessian_images = [_hessian_matrix_det(image, s) for s in sigma_list] + image_cube = np.dstack(hessian_images) + + local_maxima = peak_local_max(image_cube, threshold_abs=threshold, + footprint=np.ones((3, 3, 3)), + threshold_rel=0.0, + exclude_border=False) + + # Convert the last index to its corresponding scale value + local_maxima[:, 2] = sigma_list[local_maxima[:, 2]] + return _prune_blobs(local_maxima, overlap) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index f4894aa5..e1dd519b 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -7,6 +7,8 @@ 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 +from ._hessian_det_appx import _hessian_matrix_det +from ..transform import integral_image def _compute_derivatives(image, mode='constant', cval=0): @@ -170,6 +172,45 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0): return Hxx, Hxy, Hyy +def hessian_matrix_det(image, sigma): + """Computes the approximate Hessian Determinant over an image. + + This method uses box filters over integral images to compute the + approximate Hessian Determinant as described in [1]_. + + Parameters + ---------- + image : array + The image over which to compute Hessian Determinant. + sigma : float + Standard deviation used for the Gaussian kernel, used for the Hessian + matrix. + + Returns + ------- + out : array + The array of the Determinant of Hessians. + + References + ---------- + .. [1] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool, + "SURF: Speeded Up Robust Features" + ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf + + Notes + ----- + The running time of this method only depends on size of the image. It is + independent of `sigma` as one would expect. The downside is that the + result for `sigma` less than `3` is not accurate, i.e., not similar to + the result obtained if someone computed the Hessian and took it's + determinant. + """ + + image = img_as_float(image) + image = integral_image(image) + return np.array(_hessian_matrix_det(image, sigma)) + + 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 diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 9d8b4fbb..d7be3dcf 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None): cython(['orb_cy.pyx'], working_path=base_path) cython(['brief_cy.pyx'], working_path=base_path) cython(['_texture.pyx'], working_path=base_path) + cython(['_hessian_det_appx.pyx'], working_path=base_path) config.add_extension('corner_cy', sources=['corner_cy.c'], include_dirs=[get_numpy_include_dirs()]) @@ -28,6 +29,8 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) + config.add_extension('_hessian_det_appx', sources=['_hessian_det_appx.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py index 3ae949c2..a46c442b 100644 --- a/skimage/feature/tests/test_blob.py +++ b/skimage/feature/tests/test_blob.py @@ -1,12 +1,14 @@ import numpy as np from skimage.draw import circle -from skimage.feature import blob_dog, blob_log +from skimage.feature import blob_dog, blob_log, blob_doh import math +from numpy.testing import assert_raises def test_blob_dog(): r2 = math.sqrt(2) img = np.ones((512, 512)) + img3 = np.ones((5, 5, 5)) xs, ys = circle(400, 130, 5) img[xs, ys] = 255 @@ -18,7 +20,7 @@ def test_blob_dog(): img[xs, ys] = 255 blobs = blob_dog(img, min_sigma=5, max_sigma=50) - radius = lambda x: r2*x[2] + radius = lambda x: r2 * x[2] s = sorted(blobs, key=radius) thresh = 5 @@ -37,10 +39,13 @@ def test_blob_dog(): assert abs(b[1] - 350) <= thresh assert abs(radius(b) - 45) <= thresh + assert_raises(ValueError, blob_dog, img3) + def test_blob_log(): r2 = math.sqrt(2) img = np.ones((512, 512)) + img3 = np.ones((5, 5, 5)) xs, ys = circle(400, 130, 5) img[xs, ys] = 255 @@ -56,7 +61,7 @@ def test_blob_log(): blobs = blob_log(img, min_sigma=5, max_sigma=20, threshold=1) - radius = lambda x: r2*x[2] + radius = lambda x: r2 * x[2] s = sorted(blobs, key=radius) thresh = 3 @@ -79,3 +84,131 @@ def test_blob_log(): assert abs(b[0] - 200) <= thresh assert abs(b[1] - 350) <= thresh assert abs(radius(b) - 30) <= thresh + + # Testing log scale + blobs = blob_log( + img, + min_sigma=5, + max_sigma=20, + threshold=1, + log_scale=True) + + b = s[0] + assert abs(b[0] - 400) <= thresh + assert abs(b[1] - 130) <= thresh + assert abs(radius(b) - 5) <= thresh + + b = s[1] + assert abs(b[0] - 160) <= thresh + assert abs(b[1] - 50) <= thresh + assert abs(radius(b) - 15) <= thresh + + b = s[2] + assert abs(b[0] - 100) <= thresh + assert abs(b[1] - 300) <= thresh + assert abs(radius(b) - 25) <= thresh + + b = s[3] + assert abs(b[0] - 200) <= thresh + assert abs(b[1] - 350) <= thresh + assert abs(radius(b) - 30) <= thresh + + assert_raises(ValueError, blob_log, img3) + + +def test_blob_doh(): + img = np.ones((512, 512), dtype=np.uint8) + img3 = np.ones((5, 5, 5)) + + xs, ys = circle(400, 130, 20) + img[xs, ys] = 255 + + xs, ys = circle(460, 50, 30) + img[xs, ys] = 255 + + xs, ys = circle(100, 300, 40) + img[xs, ys] = 255 + + xs, ys = circle(200, 350, 50) + img[xs, ys] = 255 + + blobs = blob_doh( + img, + min_sigma=1, + max_sigma=60, + num_sigma=10, + threshold=.05) + + radius = lambda x: x[2] + s = sorted(blobs, key=radius) + thresh = 3 + + b = s[0] + assert abs(b[0] - 400) <= thresh + assert abs(b[1] - 130) <= thresh + assert abs(radius(b) - 20) <= thresh + + b = s[1] + assert abs(b[0] - 460) <= thresh + assert abs(b[1] - 50) <= thresh + assert abs(radius(b) - 30) <= thresh + + b = s[2] + assert abs(b[0] - 100) <= thresh + assert abs(b[1] - 300) <= thresh + assert abs(radius(b) - 40) <= thresh + + b = s[3] + assert abs(b[0] - 200) <= thresh + assert abs(b[1] - 350) <= thresh + assert abs(radius(b) - 50) <= thresh + + # Testing log scale + blobs = blob_doh( + img, + min_sigma=1, + max_sigma=60, + num_sigma=10, + log_scale=True, + threshold=.05) + + b = s[0] + assert abs(b[0] - 400) <= thresh + assert abs(b[1] - 130) <= thresh + assert abs(radius(b) - 20) <= thresh + + b = s[1] + assert abs(b[0] - 460) <= thresh + assert abs(b[1] - 50) <= thresh + assert abs(radius(b) - 30) <= thresh + + b = s[2] + assert abs(b[0] - 100) <= thresh + assert abs(b[1] - 300) <= thresh + assert abs(radius(b) - 40) <= thresh + + b = s[3] + assert abs(b[0] - 200) <= thresh + assert abs(b[1] - 350) <= thresh + assert abs(radius(b) - 50) <= thresh + + assert_raises(ValueError, blob_doh, img3) + + +def test_blob_overlap(): + img = np.ones((512, 512), dtype=np.uint8) + + xs, ys = circle(100, 100, 20) + img[xs, ys] = 255 + + xs, ys = circle(120, 100, 30) + img[xs, ys] = 255 + + blobs = blob_doh( + img, + min_sigma=1, + max_sigma=60, + num_sigma=10, + threshold=.05) + + assert len(blobs) == 1 diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index 21070779..7d2af0df 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -12,7 +12,8 @@ from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi, corner_kitchen_rosenfeld, corner_foerstner, corner_fast, corner_orientations, structure_tensor, structure_tensor_eigvals, - hessian_matrix, hessian_matrix_eigvals) + hessian_matrix, hessian_matrix_eigvals, + hessian_matrix_det) def test_structure_tensor(): @@ -91,6 +92,13 @@ def test_hessian_matrix_eigvals(): [0, 0, 0, 0, 0]])) +def test_hessian_matrix_det(): + image = np.zeros((5, 5)) + image[2, 2] = 1 + det = hessian_matrix_det(image, 5) + assert_almost_equal(det, 0, decimal = 3) + + def test_square_image(): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1.