From 7d407fc778ec0511767bf503e3d05a9824e4c9c6 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 6 Apr 2013 20:20:37 -0500 Subject: [PATCH 1/4] Fix kernel-size calculation for non-zero theta. --- skimage/filter/_gabor.py | 31 ++++++++++++++------------ skimage/filter/tests/test_gabor.py | 35 +++++++++++++++++++++++------- 2 files changed, 44 insertions(+), 22 deletions(-) diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py index 02210a8f..c224101e 100644 --- a/skimage/filter/_gabor.py +++ b/skimage/filter/_gabor.py @@ -5,16 +5,16 @@ from scipy import ndimage def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0): """Build complex 2D Gabor filter kernel. - Frequency and orientation representations of the Gabor filter are similar to - those of the human visual system. It is especially suitable for texture + Frequency and orientation representations of the Gabor filter are similar + to those of the human visual system. It is especially suitable for texture classification using Gabor filter banks. Parameters ---------- - sigma_x : float - Standard deviation in x-direction. - sigma_y : float - Standard deviation in y-direction. + 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 @@ -34,8 +34,11 @@ def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0): """ - x0 = np.ceil(max(3 * sigma_x, 1)) - y0 = np.ceil(max(3 * sigma_y, 1)) + n_stds = 3 + x0 = np.ceil(max(np.abs(n_stds * sigma_x * np.cos(theta)), + np.abs(n_stds * sigma_y * np.sin(theta)), 1)) + y0 = np.ceil(max(np.abs(n_stds * sigma_y * np.cos(theta)), + np.abs(n_stds * sigma_x * np.sin(theta)), 1)) y, x = np.mgrid[-y0:y0+1, -x0:x0+1] rotx = x * np.cos(theta) + y * np.sin(theta) @@ -56,16 +59,16 @@ def gabor_filter(image, sigma_x, sigma_y, frequency, theta, offset=0, The real and imaginary parts of the Gabor filter kernel are applied to the image. - Frequency and orientation representations of the Gabor filter are similar to - those of the human visual system. It is especially suitable for texture + Frequency and orientation representations of the Gabor filter are similar + to those of the human visual system. It is especially suitable for texture classification using Gabor filter banks. Parameters ---------- - sigma_x : float - Standard deviation in x-direction. - sigma_y : float - Standard deviation in y-direction. + 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 diff --git a/skimage/filter/tests/test_gabor.py b/skimage/filter/tests/test_gabor.py index 4080aa17..58111ce0 100644 --- a/skimage/filter/tests/test_gabor.py +++ b/skimage/filter/tests/test_gabor.py @@ -1,25 +1,44 @@ import numpy as np -from numpy.testing import assert_almost_equal, assert_array_almost_equal +from numpy.testing import (assert_equal, assert_almost_equal, + assert_array_almost_equal) from skimage.filter import gabor_kernel, gabor_filter +def test_gabor_kernel_size(): + sigma_x = 5 + sigma_y = 10 + # Sizes cut off at +/- three sigma + 1 for the center + size_x = sigma_x * 6 + 1 + size_y = sigma_y * 6 + 1 + + theta = 0 + kernel = gabor_kernel(sigma_x, sigma_y, 0, theta) + assert_equal(kernel.shape, (size_y, size_x)) + + theta = np.pi / 2 + kernel = gabor_kernel(sigma_x, sigma_y, 0, theta) + assert_equal(kernel.shape, (size_x, size_y)) + + + def test_gabor_kernel_sum(): - for sigmax in range(1, 10, 2): - for sigmay in range(1, 10, 2): + 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(sigmax, sigmay, frequency+0.1, 0) + kernel = gabor_kernel(sigma_x, sigma_y, frequency+0.1, 0) # make sure gaussian distribution is covered nearly 100% assert_almost_equal(np.abs(kernel).sum(), 1, 2) def test_gabor_kernel_theta(): - for sigmax in range(1, 10, 2): - for sigmay in range(1, 10, 2): + for sigma_x in range(1, 10, 2): + 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(sigmax, sigmay, frequency+0.1, theta) - kernel180 = gabor_kernel(sigmax, sigmay, frequency, + kernel0 = gabor_kernel(sigma_x, sigma_y, frequency+0.1, + theta) + kernel180 = gabor_kernel(sigma_x, sigma_y, frequency, theta+np.pi) assert_array_almost_equal(np.abs(kernel0), From 2bf178ceb2e59c0d9bd9d124d00cd4f5af2aa000 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 6 Apr 2013 21:18:53 -0500 Subject: [PATCH 2/4] 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__": From 67536c94f64422fa08b7735e2f3ac59def78c0b0 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 6 Apr 2013 21:30:50 -0500 Subject: [PATCH 3/4] DOC: Tweak docstrings for Gabor filter --- skimage/filter/_gabor.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py index c2ab2860..f766ac74 100644 --- a/skimage/filter/_gabor.py +++ b/skimage/filter/_gabor.py @@ -13,7 +13,7 @@ def _sigma_prefactor(bandwidth): def gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, offset=0): - """Build complex 2D Gabor filter kernel. + """Return complex 2D Gabor filter kernel. Frequency and orientation representations of the Gabor filter are similar to those of the human visual system. It is especially suitable for texture @@ -72,10 +72,10 @@ def gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, offset=0, mode='reflect', cval=0): - """Perform Gabor filtering. + """Return real and imaginary responses to Gabor filter. The real and imaginary parts of the Gabor filter kernel are applied to the - image. + image and the response is returned as a pair of arrays. Frequency and orientation representations of the Gabor filter are similar to those of the human visual system. It is especially suitable for texture @@ -83,6 +83,8 @@ def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None, Parameters ---------- + image : array + Input image. frequency : float Frequency of the harmonic function. theta : float @@ -100,7 +102,7 @@ def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None, Returns ------- - real, imag : complex arrays + real, imag : arrays Filtered images using the real and imaginary parts of the Gabor filter kernel. From ba12acdeb816a5071c26933fc02498cf05b0023b Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 6 Apr 2013 22:34:48 -0500 Subject: [PATCH 4/4] Update gabor example. The parameter order to `gabor_filter` changed so this example was broken. Also, add plots of the Gabor responses to the demo. --- doc/examples/plot_gabor.py | 65 ++++++++++++++++++++++++++------------ 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/doc/examples/plot_gabor.py b/doc/examples/plot_gabor.py index a72aa78f..b57bd353 100644 --- a/doc/examples/plot_gabor.py +++ b/doc/examples/plot_gabor.py @@ -51,22 +51,24 @@ for theta in range(4): theta = theta / 4. * np.pi for sigma in (1, 3): for frequency in (0.05, 0.25): - kernel = np.real(gabor_kernel(sigma, sigma, frequency, theta)) + kernel = np.real(gabor_kernel(frequency, theta=theta, + sigma_x=sigma, sigma_y=sigma)) kernels.append(kernel) -brick = img_as_float(data.load('brick.png')) -grass = img_as_float(data.load('grass.png')) -wall = img_as_float(data.load('rough-wall.png')) +shrink = (slice(0, None, 3), slice(0, None, 3)) +brick = img_as_float(data.load('brick.png'))[shrink] +grass = img_as_float(data.load('grass.png'))[shrink] +wall = img_as_float(data.load('rough-wall.png'))[shrink] image_names = ('brick', 'grass', 'wall') +images = (brick, grass, wall) -# prepare refernce features +# prepare reference features ref_feats = np.zeros((3, len(kernels), 2), dtype=np.double) ref_feats[0, :, :] = compute_feats(brick, kernels) ref_feats[1, :, :] = compute_feats(grass, kernels) ref_feats[2, :, :] = compute_feats(wall, kernels) - print 'Rotated images matched against references using Gabor filter banks:' print 'original: brick, rotated: 30deg, match result:', @@ -82,29 +84,50 @@ feats = compute_feats(nd.rotate(grass, angle=145, reshape=False), kernels) print image_names[match(feats, ref_feats)] -# plot a selection of the filter bank kernels +def power(image, kernel): + # Normalize images for better comparison. + image = (image - image.mean()) / image.std() + return np.sqrt(nd.convolve(image, np.real(kernel), mode='wrap')**2 + + nd.convolve(image, np.imag(kernel), mode='wrap')**2) -kernels = [] +# Plot a selection of the filter bank kernels and their responses. +results = [] kernel_params = [] -for theta in (0, 1, 3): +for theta in (0, 1): theta = theta / 4. * np.pi - for frequency in (0.05, 0.1, 0.25): - kernel = np.real(gabor_kernel(10, 10, frequency, theta)) - kernels.append(kernel) - params = 'theta=%d, frequency=%.2f' % (theta * 180 / np.pi, frequency) + for frequency in (0.1, 0.4): + kernel = gabor_kernel(frequency, theta=theta) + params = 'theta=%d,\nfrequency=%.2f' % (theta * 180 / np.pi, frequency) kernel_params.append(params) + # Save kernel and the power image for each image + results.append((kernel, [power(img, kernel) for img in images])) - -fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3, - figsize=(9, 6)) +fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(9, 6)) plt.gray() -fig.text(.5, .95, 'Gabor filter bank kernels', - horizontalalignment='center', fontsize=15) +fig.suptitle('Image responses for Gabor filter kernels', fontsize=15) -for i, ax in enumerate((ax1, ax2, ax3, ax4, ax5, ax6)): - ax.imshow(kernels[i], interpolation='nearest') +axes[0][0].axis('off') + +# Plot original images +for label, img, ax in zip(image_names, images, axes[0][1:]): + ax.imshow(img) + ax.set_title(label) ax.axis('off') - ax.set_title(kernel_params[i]) + +for label, (kernel, powers), ax_row in zip(kernel_params, results, axes[1:]): + # Plot Gabor kernel + ax = ax_row[0] + ax.imshow(np.real(kernel), interpolation='nearest') + ax.set_ylabel(label) + ax.set_xticks([]) + ax.set_yticks([]) + + # Plot Gabor responses with the contrast normalized for each filter + vmin = np.min(powers) + vmax = np.max(powers) + for patch, ax in zip(powers, ax_row[1:]): + ax.imshow(patch, vmin=vmin, vmax=vmax) + ax.axis('off') plt.show()