From 2bf178ceb2e59c0d9bd9d124d00cd4f5af2aa000 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 6 Apr 2013 21:18:53 -0500 Subject: [PATCH] Add `bandwidth` parameter to `gabor_kernel` Note that this changes the API of `gabor_kernel` and `gabor_filter`: The input parameters are rearranged and positional arguments are changed to keyword args. --- skimage/filter/_gabor.py | 46 ++++++++++++++++++++++-------- skimage/filter/tests/test_gabor.py | 36 ++++++++++++++++------- 2 files changed, 59 insertions(+), 23 deletions(-) diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py index c224101e..c2ab2860 100644 --- a/skimage/filter/_gabor.py +++ b/skimage/filter/_gabor.py @@ -2,7 +2,17 @@ import numpy as np from scipy import ndimage -def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0): +__all__ = ['gabor_kernel', 'gabor_filter'] + + +def _sigma_prefactor(bandwidth): + b = bandwidth + # See http://www.cs.rug.nl/~imaging/simplecell.html + return 1.0 / np.pi * np.sqrt(np.log(2)/2.0) * (2.0**b + 1) / (2.0**b - 1) + + +def gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, + offset=0): """Build complex 2D Gabor filter kernel. Frequency and orientation representations of the Gabor filter are similar @@ -11,14 +21,18 @@ def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0): Parameters ---------- + frequency : float + Frequency of the harmonic function. + theta : float + Orientation in radians. If 0, the harmonic is in the x-direction. + bandwidth : float + The bandwidth captured by the filter. For fixed bandwidth, `sigma_x` + and `sigma_y` will decrease with increasing frequency. This value is + ignored if `sigma_x` and `sigma_y` are set by the user. sigma_x, sigma_y : float Standard deviation in x- and y-directions. These directions apply to the kernel *before* rotation. If `theta = pi/2`, then the kernel is rotated 90 degrees so that `sigma_x` controls the *vertical* direction. - frequency : float - Frequency of the harmonic function. - theta : float - Orientation in radians. offset : float, optional Phase offset of harmonic function in radians. @@ -33,6 +47,10 @@ def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0): .. [2] http://mplab.ucsd.edu/tutorials/gabor.pdf """ + if sigma_x is None: + sigma_x = _sigma_prefactor(bandwidth) / frequency + if sigma_y is None: + sigma_y = _sigma_prefactor(bandwidth) / frequency n_stds = 3 x0 = np.ceil(max(np.abs(n_stds * sigma_x * np.cos(theta)), @@ -52,8 +70,8 @@ def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0): return g -def gabor_filter(image, sigma_x, sigma_y, frequency, theta, offset=0, - mode='reflect', cval=0): +def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None, + sigma_y=None, offset=0, mode='reflect', cval=0): """Perform Gabor filtering. The real and imaginary parts of the Gabor filter kernel are applied to the @@ -65,14 +83,18 @@ def gabor_filter(image, sigma_x, sigma_y, frequency, theta, offset=0, Parameters ---------- + frequency : float + Frequency of the harmonic function. + theta : float + Orientation in radians. If 0, the harmonic is in the x-direction. + bandwidth : float + The bandwidth captured by the filter. For fixed bandwidth, `sigma_x` + and `sigma_y` will decrease with increasing frequency. This value is + ignored if `sigma_x` and `sigma_y` are set by the user. sigma_x, sigma_y : float Standard deviation in x- and y-directions. These directions apply to the kernel *before* rotation. If `theta = pi/2`, then the kernel is rotated 90 degrees so that `sigma_x` controls the *vertical* direction. - frequency : float - Frequency of the harmonic function. - theta : float - Orientation in radians. offset : float, optional Phase offset of harmonic function in radians. @@ -89,7 +111,7 @@ def gabor_filter(image, sigma_x, sigma_y, frequency, theta, offset=0, """ - g = gabor_kernel(sigma_x, sigma_y, frequency, theta, offset) + g = gabor_kernel(frequency, theta, bandwidth, sigma_x, sigma_y, offset) filtered_real = ndimage.convolve(image, np.real(g), mode=mode, cval=cval) filtered_imag = ndimage.convolve(image, np.imag(g), mode=mode, cval=cval) diff --git a/skimage/filter/tests/test_gabor.py b/skimage/filter/tests/test_gabor.py index 58111ce0..68a9c967 100644 --- a/skimage/filter/tests/test_gabor.py +++ b/skimage/filter/tests/test_gabor.py @@ -2,7 +2,7 @@ import numpy as np from numpy.testing import (assert_equal, assert_almost_equal, assert_array_almost_equal) -from skimage.filter import gabor_kernel, gabor_filter +from skimage.filter._gabor import gabor_kernel, gabor_filter, _sigma_prefactor def test_gabor_kernel_size(): @@ -12,21 +12,35 @@ def test_gabor_kernel_size(): size_x = sigma_x * 6 + 1 size_y = sigma_y * 6 + 1 - theta = 0 - kernel = gabor_kernel(sigma_x, sigma_y, 0, theta) + kernel = gabor_kernel(0, theta=0, sigma_x=sigma_x, sigma_y=sigma_y) assert_equal(kernel.shape, (size_y, size_x)) - theta = np.pi / 2 - kernel = gabor_kernel(sigma_x, sigma_y, 0, theta) + kernel = gabor_kernel(0, theta=np.pi/2, sigma_x=sigma_x, sigma_y=sigma_y) assert_equal(kernel.shape, (size_x, size_y)) +def test_gabor_kernel_bandwidth(): + kernel = gabor_kernel(1, bandwidth=1) + assert_equal(kernel.shape, (5, 5)) + + kernel = gabor_kernel(1, bandwidth=0.5) + assert_equal(kernel.shape, (9, 9)) + + kernel = gabor_kernel(0.5, bandwidth=1) + assert_equal(kernel.shape, (9, 9)) + + +def test_sigma_prefactor(): + assert_almost_equal(_sigma_prefactor(1), 0.56, 2) + assert_almost_equal(_sigma_prefactor(0.5), 1.09, 2) + def test_gabor_kernel_sum(): for sigma_x in range(1, 10, 2): for sigma_y in range(1, 10, 2): for frequency in range(0, 10, 2): - kernel = gabor_kernel(sigma_x, sigma_y, frequency+0.1, 0) + kernel = gabor_kernel(frequency+0.1, theta=0, + sigma_x=sigma_x, sigma_y=sigma_y) # make sure gaussian distribution is covered nearly 100% assert_almost_equal(np.abs(kernel).sum(), 1, 2) @@ -36,17 +50,17 @@ def test_gabor_kernel_theta(): for sigma_y in range(1, 10, 2): for frequency in range(0, 10, 2): for theta in range(0, 10, 2): - kernel0 = gabor_kernel(sigma_x, sigma_y, frequency+0.1, - theta) - kernel180 = gabor_kernel(sigma_x, sigma_y, frequency, - theta+np.pi) + kernel0 = gabor_kernel(frequency+0.1, theta=theta, + sigma_x=sigma_x, sigma_y=sigma_y) + kernel180 = gabor_kernel(frequency, theta=theta+np.pi, + sigma_x=sigma_x, sigma_y=sigma_y) assert_array_almost_equal(np.abs(kernel0), np.abs(kernel180)) def test_gabor_filter(): - real, imag = gabor_filter(np.random.random((100, 100)), 1, 1, 1, 1) + real, imag = gabor_filter(np.random.random((100, 100)), 1) if __name__ == "__main__":