diff --git a/doc/examples/plot_gabor.py b/doc/examples/plot_gabor.py new file mode 100644 index 00000000..b57bd353 --- /dev/null +++ b/doc/examples/plot_gabor.py @@ -0,0 +1,133 @@ +""" +============================================= +Gabor filter banks for texture classification +============================================= + +In this example, we will see how to classify textures based on Gabor filter +banks. Frequency and orientation representations of the Gabor filter are similar +to those of the human visual system. + +The images are filtered using the real parts of various different Gabor filter +kernels. The mean and variance of the filtered images are then used as features +for classification, which is based on the least squared error for simplicity. + +""" + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from scipy import ndimage as nd +from skimage import data +from skimage.util import img_as_float +from skimage.filter import gabor_kernel + + +matplotlib.rcParams['font.size'] = 9 + + +def compute_feats(image, kernels): + feats = np.zeros((len(kernels), 2), dtype=np.double) + for k, kernel in enumerate(kernels): + filtered = nd.convolve(image, kernel, mode='wrap') + feats[k, 0] = filtered.mean() + feats[k, 1] = filtered.var() + return feats + + +def match(feats, ref_feats): + min_error = np.inf + min_i = None + for i in range(ref_feats.shape[0]): + error = np.sum((feats - ref_feats[i, :])**2) + if error < min_error: + min_error = error + min_i = i + return min_i + + +# prepare filter bank kernels +kernels = [] +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(frequency, theta=theta, + sigma_x=sigma, sigma_y=sigma)) + kernels.append(kernel) + + +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 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:', +feats = compute_feats(nd.rotate(brick, angle=190, reshape=False), kernels) +print image_names[match(feats, ref_feats)] + +print 'original: brick, rotated: 70deg, match result:', +feats = compute_feats(nd.rotate(brick, angle=70, reshape=False), kernels) +print image_names[match(feats, ref_feats)] + +print 'original: grass, rotated: 145deg, match result:', +feats = compute_feats(nd.rotate(grass, angle=145, reshape=False), kernels) +print image_names[match(feats, ref_feats)] + + +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) + +# Plot a selection of the filter bank kernels and their responses. +results = [] +kernel_params = [] +for theta in (0, 1): + theta = theta / 4. * np.pi + 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, axes = plt.subplots(nrows=5, ncols=4, figsize=(9, 6)) +plt.gray() + +fig.suptitle('Image responses for Gabor filter kernels', fontsize=15) + +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') + +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() diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index f32c3b23..15eec3da 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -6,4 +6,5 @@ from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt, from ._denoise import denoise_tv_chambolle, tv_denoise from ._denoise_cy import denoise_bilateral, denoise_tv_bregman from ._rank_order import rank_order +from ._gabor import gabor_kernel, gabor_filter from .thresholding import threshold_otsu, threshold_adaptive diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py new file mode 100644 index 00000000..f766ac74 --- /dev/null +++ b/skimage/filter/_gabor.py @@ -0,0 +1,121 @@ +import numpy as np +from scipy import ndimage + + +__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): + """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 + classification using Gabor filter banks. + + 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. + offset : float, optional + Phase offset of harmonic function in radians. + + Returns + ------- + g : complex array + Complex filter kernel. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Gabor_filter + .. [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)), + 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) + roty = -x * np.sin(theta) + y * np.cos(theta) + + g = np.zeros(y.shape, dtype=np.complex) + g[:] = np.exp(-0.5 * (rotx**2 / sigma_x**2 + roty**2 / sigma_y**2)) + g /= 2 * np.pi * sigma_x * sigma_y + g *= np.exp(1j * (2 * np.pi * frequency * rotx + offset)) + + return g + + +def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None, + sigma_y=None, offset=0, mode='reflect', cval=0): + """Return real and imaginary responses to Gabor filter. + + The real and imaginary parts of the Gabor filter kernel are applied to the + 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 + classification using Gabor filter banks. + + Parameters + ---------- + image : array + Input image. + 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. + offset : float, optional + Phase offset of harmonic function in radians. + + Returns + ------- + real, imag : arrays + Filtered images using the real and imaginary parts of the Gabor filter + kernel. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Gabor_filter + .. [2] http://mplab.ucsd.edu/tutorials/gabor.pdf + + """ + + 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) + + return filtered_real, filtered_imag diff --git a/skimage/filter/tests/test_gabor.py b/skimage/filter/tests/test_gabor.py new file mode 100644 index 00000000..7035d82a --- /dev/null +++ b/skimage/filter/tests/test_gabor.py @@ -0,0 +1,82 @@ +import numpy as np +from numpy.testing import (assert_equal, assert_almost_equal, + assert_array_almost_equal) + +from skimage.filter._gabor import gabor_kernel, gabor_filter, _sigma_prefactor + + +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 + + kernel = gabor_kernel(0, theta=0, sigma_x=sigma_x, sigma_y=sigma_y) + assert_equal(kernel.shape, (size_y, size_x)) + + 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(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) + + +def test_gabor_kernel_theta(): + 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(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(): + Y, X = np.mgrid[:40, :40] + frequencies = (0.1, 0.3) + wave_images = [np.sin(2 * np.pi * X * f) for f in frequencies] + + def match_score(image, frequency): + gabor_responses = gabor_filter(image, frequency) + return np.mean(np.hypot(*gabor_responses)) + + # Gabor scores: diagonals are frequency-matched, off-diagonals are not. + responses = np.array([[match_score(image, f) for f in frequencies] + for image in wave_images]) + assert responses[0, 0] > responses[0, 1] + assert responses[1, 1] > responses[0, 1] + assert responses[0, 0] > responses[1, 0] + assert responses[1, 1] > responses[1, 0] + + +if __name__ == "__main__": + from numpy import testing + testing.run_module_suite()