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

Merge pull request #930 from vighneshbirodkar/doh

Browse Source 
Blob Detection - Determinant of Hessian
This commit is contained in:
Johannes L. Schönberger
2014-04-11 07:33:16 -04:00
parent 88358f67cf 4d11868c76
commit 0c13cce2dd
8 changed files with 463 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files
bento.info
+3
View File
@@ -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
skimage/feature/__init__.py
+4 -2
View File
@@ -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']
skimage/feature/_hessian_det_appx.pyx
+153
View File
@@ -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
skimage/feature/blob.py
+114 -4
View File
@@ -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)
skimage/feature/corner.py
+41
View File
@@ -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
skimage/feature/setup.py
+3
View File
@@ -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
skimage/feature/tests/test_blob.py
+136 -3
View File
@@ -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
skimage/feature/tests/test_corner.py
+9 -1
View File
@@ -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.