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

Merge pull request #2 from tonysyu/gabor

Browse Source 
Change API for Gabor filter
This commit is contained in:
Johannes Schönberger
2013-04-07 00:22:58 -07:00
parent 95d1e627c6 ba12acdeb8
commit 5c45812d92
3 changed files with 139 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_gabor.py
+44 -21
View File
@@ -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()
skimage/filter/_gabor.py
+51 -24
View File
@@ -2,23 +2,37 @@ import numpy as np
from scipy import ndimage
def gabor_kernel(sigma_x, sigma_y, frequency, theta, offset=0):
"""Build complex 2D Gabor filter kernel.
__all__ = ['gabor_kernel', 'gabor_filter']
Frequency and orientation representations of the Gabor filter are similar to
those of the human visual system. It is especially suitable for texture
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
----------
sigma_x : float
Standard deviation in x-direction.
sigma_y : float
Standard deviation in y-direction.
frequency : float
Frequency of the harmonic function.
theta : float
Orientation in radians.
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.
@@ -33,9 +47,16 @@ 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
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)
@@ -49,33 +70,39 @@ 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):
"""Perform Gabor filtering.
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.
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
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.
image : array
Input image.
frequency : float
Frequency of the harmonic function.
theta : float
Orientation in radians.
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 : complex arrays
real, imag : arrays
Filtered images using the real and imaginary parts of the Gabor filter
kernel.
@@ -86,7 +113,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)
skimage/filter/tests/test_gabor.py
+44 -11
View File
@@ -1,33 +1,66 @@
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
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 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(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 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,
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__":