mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-13 12:40:24 +08:00
Tony S Yu
@@ -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()
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user