diff --git a/doc/examples/filters/plot_denoise.py b/doc/examples/filters/plot_denoise.py
index e66d2478..be9314ee 100644
--- a/doc/examples/filters/plot_denoise.py
+++ b/doc/examples/filters/plot_denoise.py
@@ -49,14 +49,14 @@ ax[0, 0].set_title('noisy')
 ax[0, 1].imshow(denoise_tv_chambolle(noisy, weight=0.1, multichannel=True))
 ax[0, 1].axis('off')
 ax[0, 1].set_title('TV')
-ax[0, 2].imshow(denoise_bilateral(noisy, sigma_range=0.05, sigma_spatial=15))
+ax[0, 2].imshow(denoise_bilateral(noisy, sigma_range=0.05, sigma_spatial=15, multichannel=True))
 ax[0, 2].axis('off')
 ax[0, 2].set_title('Bilateral')
 
 ax[1, 0].imshow(denoise_tv_chambolle(noisy, weight=0.2, multichannel=True))
 ax[1, 0].axis('off')
 ax[1, 0].set_title('(more) TV')
-ax[1, 1].imshow(denoise_bilateral(noisy, sigma_range=0.1, sigma_spatial=15))
+ax[1, 1].imshow(denoise_bilateral(noisy, sigma_range=0.1, sigma_spatial=15, multichannel=True))
 ax[1, 1].axis('off')
 ax[1, 1].set_title('(more) Bilateral')
 ax[1, 2].imshow(astro)
diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 754fcf97..7c172744 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -46,7 +46,7 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
     cval : string
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
-    multichannel : bool, default False
+    multichannel : bool
         Whether the last axis of the image is to be interpreted as multiple
         channels or another spatial dimension.
 
@@ -66,17 +66,27 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
     >>> astro = astro[220:300, 220:320]
     >>> noisy = astro + 0.6 * astro.std() * np.random.random(astro.shape)
     >>> noisy = np.clip(noisy, 0, 1)
-    >>> denoised = denoise_bilateral(noisy, sigma_range=0.05, sigma_spatial=15)
+    >>> denoised = denoise_bilateral(noisy, sigma_range=0.05,
+    ...                              sigma_spatial=15, multichannel=True)
     """
     if multichannel:
-        if image.shape[2] not in (3, 4):
-            msg = "Input image must be grayscale, RGB, or RGBA; but has " \
-                  "a shape {0}."
-            warnings.warn(msg.format(image.shape))
+        if image.ndim != 3:
+            raise ValueError("Use ``multichannel=False`` for 2D grayscale images. "
+                             "The last axis of the input image must be multiple "
+                             "color channels not another spatial dimension.")
+        elif image.shape[2] not in (3, 4):
+            if image.shape[2] > 4:
+                warnings.warn("The last axis of the input image is interpreted "
+                              "as channels. Input image with shape {0} has {1} "
+                              "channels in last axis. ``denoise_bilateral`` is "
+                              "implemented for 2D grayscale and color images "
+                              "only.".format(image.shape, image.shape[2]))
+            else:
+                msg = "Input image must be grayscale, RGB, or RGBA; but has shape {0}."
+                warnings.warn(msg.format(image.shape))
     else:
         if image.ndim > 2:
-            msg = "Input image must be grayscale, RGB, or RGBA; but has " \
-                  "a shape {0}."
+            msg = "Input image must be grayscale, RGB, or RGBA; but has shape {0}."
             raise TypeError(msg.format(image.shape))
 
 
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 58f0a261..8f9b8649 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -175,15 +175,19 @@ def test_denoise_bilateral_3d():
     img = np.clip(img, 0, 1)
 
     out1 = restoration.denoise_bilateral(img, sigma_range=0.1,
-                                         sigma_spatial=20)
+                                         sigma_spatial=20, multichannel=True)
     out2 = restoration.denoise_bilateral(img, sigma_range=0.2,
-                                         sigma_spatial=30)
+                                         sigma_spatial=30, multichannel=True)
 
     # make sure noise is reduced in the checkerboard cells
     assert img[30:45, 5:15].std() > out1[30:45, 5:15].std()
     assert out1[30:45, 5:15].std() > out2[30:45, 5:15].std()
 
 
+def test_denoise_bilateral_3d_grayscale():
+    img =  np.ones((500, 500, 3))
+    assert_raises(TypeError, restoration.denoise_bilateral, img)
+
 def test_denoise_bilateral_nan():
     img = np.NaN + np.empty((50, 50))
     out = restoration.denoise_bilateral(img)