From 5f70bbd5616c62a97ab75e54c17f7a5250ffaff1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 28 Oct 2012 20:00:10 +0100
Subject: [PATCH 01/11] Add gabor filter function

---
 skimage/filter/__init__.py |  1 +
 skimage/filter/_gabor.py   | 94 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 95 insertions(+)
 create mode 100644 skimage/filter/_gabor.py

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index f1c1fd49..2a957a62 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 tv_denoise, denoise_tv
 from ._denoise import denoise_bilateral
 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..8c2e107a
--- /dev/null
+++ b/skimage/filter/_gabor.py
@@ -0,0 +1,94 @@
+import numpy as np
+from scipy import ndimage
+
+
+def gabor_kernel(sigmax, sigmay, 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
+    classification using Gabor filter banks.
+
+    Parameters
+    ----------
+    sigmax : float
+        Standard deviation in x-direction.
+    sigmay : float
+        Standard deviation in y-direction.
+    frequency : float
+        Frequency of the harmonic function.
+    theta : float
+        Orientation in radians.
+    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
+
+    """
+
+    x0 = np.ceil(max(3 * sigmax, 1))
+    y0 = np.ceil(max(3 * sigmay, 1))
+    y, x = np.mgrid[-x0:x0+1, -y0:y0+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 / sigmax**2 + roty**2 / sigmay**2))
+    g /= 2 * np.pi * sigmax * sigmay
+    g *= np.exp(1j * (2 * np.pi * frequency * rotx + offset))
+
+    return g
+
+
+def gabor_filter(image, sigmax, sigmay, frequency, theta, offset=0,
+                 mode='reflect', cval=0):
+    """Perform Gabor filtering.
+
+    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
+    classification using Gabor filter banks.
+
+    Parameters
+    ----------
+    sigmax : float
+        Standard deviation in x-direction.
+    sigmay : float
+        Standard deviation in y-direction.
+    frequency : float
+        Frequency of the harmonic function.
+    theta : float
+        Orientation in radians.
+    offset : float, optional
+        Phase offset of harmonic function in radians.
+
+    Returns
+    -------
+    real, imag : complex 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(sigmax, sigmay, frequency, theta, offset=0)
+
+    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

From d240e9141365779b86e10a1aeed52bac903b9fb8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 28 Oct 2012 22:32:09 +0100
Subject: [PATCH 02/11] Fix meshgrid bug in gabor kernel and argument passing

---
 skimage/filter/_gabor.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py
index 8c2e107a..3424e12f 100644
--- a/skimage/filter/_gabor.py
+++ b/skimage/filter/_gabor.py
@@ -36,7 +36,7 @@ def gabor_kernel(sigmax, sigmay, frequency, theta, offset=0):
 
     x0 = np.ceil(max(3 * sigmax, 1))
     y0 = np.ceil(max(3 * sigmay, 1))
-    y, x = np.mgrid[-x0:x0+1, -y0:y0+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)
@@ -86,7 +86,7 @@ def gabor_filter(image, sigmax, sigmay, frequency, theta, offset=0,
 
     """
 
-    g = gabor_kernel(sigmax, sigmay, frequency, theta, offset=0)
+    g = gabor_kernel(sigmax, sigmay, frequency, theta, 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)

From aeaf2e1a1207f2094ea4298b1ecff015f5996b5a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 28 Oct 2012 22:32:28 +0100
Subject: [PATCH 03/11] Add test cases for gabor filter

---
 skimage/filter/tests/test_gabor.py | 35 ++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)
 create mode 100644 skimage/filter/tests/test_gabor.py

diff --git a/skimage/filter/tests/test_gabor.py b/skimage/filter/tests/test_gabor.py
new file mode 100644
index 00000000..4080aa17
--- /dev/null
+++ b/skimage/filter/tests/test_gabor.py
@@ -0,0 +1,35 @@
+import numpy as np
+from numpy.testing import assert_almost_equal, assert_array_almost_equal
+
+from skimage.filter import gabor_kernel, gabor_filter
+
+
+def test_gabor_kernel_sum():
+    for sigmax in range(1, 10, 2):
+        for sigmay in range(1, 10, 2):
+            for frequency in range(0, 10, 2):
+                kernel = gabor_kernel(sigmax, sigmay, 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 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)
+
+                    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)
+
+
+if __name__ == "__main__":
+    from numpy import testing
+    testing.run_module_suite()

From f5795f15bed0d90aed885e00a3b59c5ee1e0be2f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Mon, 29 Oct 2012 00:16:22 +0100
Subject: [PATCH 04/11] Add example script for gabor filter

---
 doc/examples/plot_gabor.py | 114 +++++++++++++++++++++++++++++++++++++
 1 file changed, 114 insertions(+)
 create mode 100644 doc/examples/plot_gabor.py

diff --git a/doc/examples/plot_gabor.py b/doc/examples/plot_gabor.py
new file mode 100644
index 00000000..624be9f2
--- /dev/null
+++ b/doc/examples/plot_gabor.py
@@ -0,0 +1,114 @@
+"""
+=============================================
+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 = []
+kernel_params = []
+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))
+            kernels.append(kernel)
+            params = 'theta=%d, sigma=%d,\nfrequency=%.2f' % (
+                      theta * 180 / np.pi, sigma, frequency)
+            kernel_params.append(params)
+
+
+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'))
+image_names = ('brick', 'grass', 'wall')
+
+# prepare refernce 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)]
+
+
+# plot a selection of the filter bank kernels
+
+kernels = []
+kernel_params = []
+for theta in (0, 1, 3):
+    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)
+        kernel_params.append(params)
+
+
+fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3,
+                                                       figsize=(9, 6))
+plt.gray()
+
+fig.text(.5, .95, 'Gabor filter bank kernels',
+         horizontalalignment='center', fontsize=15)
+
+for i, ax in enumerate((ax1, ax2, ax3, ax4, ax5, ax6)):
+    ax.imshow(kernels[i], interpolation='nearest')
+    ax.axis('off')
+    ax.set_title(kernel_params[i])
+
+plt.show()

From e30fa821d97ed8098ebb094f3ac80aa0d2f46bfb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Mon, 29 Oct 2012 17:46:15 +0100
Subject: [PATCH 05/11] Remove unnecessary kernel description for filtering

---
 doc/examples/plot_gabor.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/doc/examples/plot_gabor.py b/doc/examples/plot_gabor.py
index 624be9f2..a72aa78f 100644
--- a/doc/examples/plot_gabor.py
+++ b/doc/examples/plot_gabor.py
@@ -47,16 +47,12 @@ def match(feats, ref_feats):
 
 # prepare filter bank kernels
 kernels = []
-kernel_params = []
 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))
             kernels.append(kernel)
-            params = 'theta=%d, sigma=%d,\nfrequency=%.2f' % (
-                      theta * 180 / np.pi, sigma, frequency)
-            kernel_params.append(params)
 
 
 brick = img_as_float(data.load('brick.png'))

From 95d1e627c61c21975a0c93980b420dfc484f779c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jschoenberger@demuc.de>
Date: Sat, 6 Apr 2013 19:21:08 +0200
Subject: [PATCH 06/11] Rename sigma parameters by adding an underscore as
 separator

---
 skimage/filter/_gabor.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py
index 3424e12f..02210a8f 100644
--- a/skimage/filter/_gabor.py
+++ b/skimage/filter/_gabor.py
@@ -2,7 +2,7 @@ import numpy as np
 from scipy import ndimage
 
 
-def gabor_kernel(sigmax, sigmay, frequency, theta, offset=0):
+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
@@ -11,9 +11,9 @@ def gabor_kernel(sigmax, sigmay, frequency, theta, offset=0):
 
     Parameters
     ----------
-    sigmax : float
+    sigma_x : float
         Standard deviation in x-direction.
-    sigmay : float
+    sigma_y : float
         Standard deviation in y-direction.
     frequency : float
         Frequency of the harmonic function.
@@ -34,22 +34,22 @@ def gabor_kernel(sigmax, sigmay, frequency, theta, offset=0):
 
     """
 
-    x0 = np.ceil(max(3 * sigmax, 1))
-    y0 = np.ceil(max(3 * sigmay, 1))
+    x0 = np.ceil(max(3 * sigma_x, 1))
+    y0 = np.ceil(max(3 * sigma_y, 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 / sigmax**2 + roty**2 / sigmay**2))
-    g /= 2 * np.pi * sigmax * sigmay
+    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, sigmax, sigmay, frequency, theta, offset=0,
+def gabor_filter(image, sigma_x, sigma_y, frequency, theta, offset=0,
                  mode='reflect', cval=0):
     """Perform Gabor filtering.
 
@@ -62,9 +62,9 @@ def gabor_filter(image, sigmax, sigmay, frequency, theta, offset=0,
 
     Parameters
     ----------
-    sigmax : float
+    sigma_x : float
         Standard deviation in x-direction.
-    sigmay : float
+    sigma_y : float
         Standard deviation in y-direction.
     frequency : float
         Frequency of the harmonic function.
@@ -86,7 +86,7 @@ def gabor_filter(image, sigmax, sigmay, frequency, theta, offset=0,
 
     """
 
-    g = gabor_kernel(sigmax, sigmay, frequency, theta, offset)
+    g = gabor_kernel(sigma_x, sigma_y, frequency, theta, 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)

From 7d407fc778ec0511767bf503e3d05a9824e4c9c6 Mon Sep 17 00:00:00 2001
From: Tony S Yu <tsyu80@gmail.com>
Date: Sat, 6 Apr 2013 20:20:37 -0500
Subject: [PATCH 07/11] 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 <tsyu80@gmail.com>
Date: Sat, 6 Apr 2013 21:18:53 -0500
Subject: [PATCH 08/11] 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 <tsyu80@gmail.com>
Date: Sat, 6 Apr 2013 21:30:50 -0500
Subject: [PATCH 09/11] 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 <tsyu80@gmail.com>
Date: Sat, 6 Apr 2013 22:34:48 -0500
Subject: [PATCH 10/11] 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()

From 0f2628c7db9f7102f112ac72139c21b1085195d9 Mon Sep 17 00:00:00 2001
From: Tony S Yu <tsyu80@gmail.com>
Date: Sun, 7 Apr 2013 10:43:43 -0500
Subject: [PATCH 11/11] Add better test of gabor filter

---
 skimage/filter/tests/test_gabor.py | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/skimage/filter/tests/test_gabor.py b/skimage/filter/tests/test_gabor.py
index 68a9c967..7035d82a 100644
--- a/skimage/filter/tests/test_gabor.py
+++ b/skimage/filter/tests/test_gabor.py
@@ -60,7 +60,21 @@ def test_gabor_kernel_theta():
 
 
 def test_gabor_filter():
-    real, imag = gabor_filter(np.random.random((100, 100)), 1)
+    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__":