diff --git a/doc/examples/filters/plot_denoise.py b/doc/examples/filters/plot_denoise.py
index be9314ee..9d3d0485 100644
--- a/doc/examples/filters/plot_denoise.py
+++ b/doc/examples/filters/plot_denoise.py
@@ -49,14 +49,16 @@ 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, multichannel=True))
+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, multichannel=True))
+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 7c172744..b21501c1 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -7,7 +7,7 @@ import warnings
 
 
 def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
-                      bins=10000, mode='constant', cval=0, multichannel=False):
+                      bins=10000, mode='constant', cval=0, multichannel=True):
     """Denoise image using bilateral filter.
 
     This is an edge-preserving and noise reducing denoising filter. It averages
@@ -86,8 +86,11 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
                 warnings.warn(msg.format(image.shape))
     else:
         if image.ndim > 2:
-            msg = "Input image must be grayscale, RGB, or RGBA; but has shape {0}."
-            raise TypeError(msg.format(image.shape))
+            raise TypeError("Bilateral filter is not implemented for "
+                            "grayscale images of 3 or more dimensions, "
+                            "but input image has {0} dimension. Use "
+                            "``multichannel=True`` for 2-D RGB "
+                            "images.".format(image.shape))
 
 
     mode = _mode_deprecations(mode)
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 8f9b8649..7c6f0c57 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -2,6 +2,7 @@ import numpy as np
 from numpy.testing import run_module_suite, assert_raises, assert_equal
 
 from skimage import restoration, data, color, img_as_float, measure
+from skimage._shared._warnings import expected_warnings
 
 
 np.random.seed(1234)
@@ -159,16 +160,16 @@ def test_denoise_bilateral_2d():
     img = np.clip(img, 0, 1)
 
     out1 = restoration.denoise_bilateral(img, sigma_range=0.1,
-                                         sigma_spatial=20)
+                                         sigma_spatial=20, multichannel=False)
     out2 = restoration.denoise_bilateral(img, sigma_range=0.2,
-                                         sigma_spatial=30)
+                                         sigma_spatial=30, multichannel=False)
 
     # 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():
+def test_denoise_bilateral_color():
     img = checkerboard.copy()
     # add some random noise
     img += 0.5 * img.std() * np.random.rand(*img.shape)
@@ -185,12 +186,31 @@ def test_denoise_bilateral_3d():
 
 
 def test_denoise_bilateral_3d_grayscale():
-    img =  np.ones((500, 500, 3))
-    assert_raises(TypeError, restoration.denoise_bilateral, img)
+    img =  np.ones((50, 50, 3))
+    assert_raises(TypeError, restoration.denoise_bilateral, img, \
+                  multichannel=False)
+
+
+def test_denoise_bilateral_3d_multichannel():
+    img = np.ones((50, 50, 50))
+    with expected_warnings(["grayscale"]):
+        result = restoration.denoise_bilateral(img, multichannel=True)
+
+    expected = np.empty_like(img)
+    expected.fill(np.nan)
+
+    assert_equal(result, expected)
+
+
+def test_denoise_bilateral_multidimensional():
+    img = np.ones((10, 10, 10, 10))
+    assert_raises(ValueError, restoration.denoise_bilateral, img,
+                  multichannel=True)
+
 
 def test_denoise_bilateral_nan():
     img = np.NaN + np.empty((50, 50))
-    out = restoration.denoise_bilateral(img)
+    out = restoration.denoise_bilateral(img, multichannel=False)
     assert_equal(img, out)