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),