diff --git a/skimage/_shared/utils.py b/skimage/_shared/utils.py
index da64a32d..49cb3294 100644
--- a/skimage/_shared/utils.py
+++ b/skimage/_shared/utils.py
@@ -143,8 +143,23 @@ def safe_as_int(val, atol=1e-3):
     return np.round(val).astype(np.int64)
 
 
-def assert_nD(array, arg_name='image', ndim=2):
+def assert_nD(array, ndim, arg_name='image'):
+    """
+    Verify an array meets the desired ndims.
+
+    Parameters
+    ----------
+    array : array-like
+        Input array to be validated
+    ndim : int or iterable of ints
+        Allowable ndim or ndims for the array.
+    arg_name : str, optional
+        The name of the array in the original function.
+
+    """
     array = np.asanyarray(array)
-    if array.ndim != ndim:
-        msg = "The parameter `%s` must be a %s-dimensional array"
-        raise ValueError(msg % (arg_name, ndim))
+    msg = "The parameter `%s` must be a %s-dimensional array"
+    if isinstance(ndim, int):
+        ndim = [ndim]
+    if not array.ndim in ndim:
+        raise ValueError(msg % (arg_name, '-or-'.join([str(n) for n in ndim])))
diff --git a/skimage/feature/_canny.py b/skimage/feature/_canny.py
index 8dcd1570..a1ba2fb8 100644
--- a/skimage/feature/_canny.py
+++ b/skimage/feature/_canny.py
@@ -17,6 +17,7 @@ import scipy.ndimage as ndi
 from scipy.ndimage import (gaussian_filter,
                            generate_binary_structure, binary_erosion, label)
 from skimage import dtype_limits
+from skimage._shared.utils import assert_nD
 
 
 def smooth_with_function_and_mask(image, function, mask):
@@ -148,9 +149,7 @@ def canny(image, sigma=1., low_threshold=None, high_threshold=None, mask=None):
     # mask by one and then mask the output. We also mask out the border points
     # because who knows what lies beyond the edge of the image?
     #
-
-    if image.ndim != 2:
-        raise TypeError("The input 'image' must be a two-dimensional array.")
+    assert_nD(image, 2)
 
     if low_threshold is None:
         low_threshold = 0.1 * dtype_limits(image)[1]
diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
index 3a55faf1..d2d03ec7 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -3,6 +3,7 @@ from scipy import sqrt, pi, arctan2, cos, sin, exp
 from scipy.ndimage import gaussian_filter
 import skimage.color
 from skimage import img_as_float, draw
+from skimage._shared.utils import assert_nD
 
 
 def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
@@ -93,9 +94,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
     .. [2] http://cvlab.epfl.ch/alumni/tola/daisy.html
     '''
 
-    # Validate image format.
-    if img.ndim != 2:
-        raise ValueError('Only grey-level images are supported.')
+    assert_nD(img, 2, 'img')
 
     img = img_as_float(img)
 
diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py
index 72411308..dd532b2e 100644
--- a/skimage/feature/_hog.py
+++ b/skimage/feature/_hog.py
@@ -1,6 +1,7 @@
 import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin
 from scipy.ndimage import uniform_filter
+from skimage._shared.utils import assert_nD
 
 
 def hog(image, orientations=9, pixels_per_cell=(8, 8),
@@ -59,8 +60,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
     shadowing and illumination variations.
     """
 
-    if image.ndim > 2:
-        raise ValueError("Currently only supports grey-level images")
+    assert_nD(image, 2)
 
     if normalise:
         image = sqrt(image)
@@ -79,7 +79,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
         # convert uint image to float
         # to avoid problems with subtracting unsigned numbers in np.diff()
         image = image.astype('float')
-    
+
     gx = np.empty(image.shape, dtype=np.double)
     gx[:, 0] = 0
     gx[:, -1] = 0
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 134e7026..6431ed90 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -9,6 +9,7 @@ from skimage.util import img_as_float
 from .peak import peak_local_max
 from ._hessian_det_appx import _hessian_matrix_det
 from skimage.transform import integral_image
+from skimage._shared.utils import assert_nD
 
 
 # This basic blob detection algorithm is based on:
@@ -169,9 +170,7 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=2.0,
     -----
     The radius of each blob is approximately :math:`\sqrt{2}sigma`.
     """
-
-    if image.ndim != 2:
-        raise ValueError("'image' must be a grayscale ")
+    assert_nD(image, 2)
 
     image = img_as_float(image)
 
@@ -275,8 +274,7 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
     The radius of each blob is approximately :math:`\sqrt{2}sigma`.
     """
 
-    if image.ndim != 2:
-        raise ValueError("'image' must be a grayscale ")
+    assert_nD(image, 2)
 
     image = img_as_float(image)
 
@@ -385,8 +383,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,
     due to the box filters used in the approximation of Hessian Determinant.
     """
 
-    if image.ndim != 2:
-        raise ValueError("'image' must be grayscale ")
+    assert_nD(image, 2)
 
     image = img_as_float(image)
     image = integral_image(image)
diff --git a/skimage/feature/brief.py b/skimage/feature/brief.py
index d1626f17..856fd6b7 100644
--- a/skimage/feature/brief.py
+++ b/skimage/feature/brief.py
@@ -5,6 +5,7 @@ from .util import (DescriptorExtractor, _mask_border_keypoints,
                    _prepare_grayscale_input_2D)
 
 from .brief_cy import _brief_loop
+from skimage._shared.utils import assert_nD
 
 
 class BRIEF(DescriptorExtractor):
@@ -137,6 +138,7 @@ class BRIEF(DescriptorExtractor):
             Keypoint coordinates as ``(row, col)``.
 
         """
+        assert_nD(image, 2)
 
         np.random.seed(self.sample_seed)
 
diff --git a/skimage/feature/censure.py b/skimage/feature/censure.py
index eb69f115..1b1754a7 100644
--- a/skimage/feature/censure.py
+++ b/skimage/feature/censure.py
@@ -9,7 +9,7 @@ from skimage.morphology import octagon, star
 from skimage.feature.util import _mask_border_keypoints
 
 from skimage.feature.censure_cy import _censure_dob_loop
-
+from skimage._shared.utils import assert_nD
 
 # The paper(Reference [1]) mentions the sizes of the Octagon shaped filter
 # kernel for the first seven scales only. The sizes of the later scales
@@ -231,6 +231,8 @@ class CENSURE(FeatureDetector):
         # (4) Finally, we remove the border keypoints and return the keypoints
         # along with its corresponding scale.
 
+        assert_nD(image, 2)
+
         num_scales = self.max_scale - self.min_scale
 
         image = np.ascontiguousarray(_prepare_grayscale_input_2D(image))
diff --git a/skimage/feature/orb.py b/skimage/feature/orb.py
index 2ddcf4f3..0ccab314 100644
--- a/skimage/feature/orb.py
+++ b/skimage/feature/orb.py
@@ -7,6 +7,7 @@ from skimage.feature.util import (FeatureDetector, DescriptorExtractor,
 from skimage.feature import (corner_fast, corner_orientations, corner_peaks,
                              corner_harris)
 from skimage.transform import pyramid_gaussian
+from skimage._shared.utils import assert_nD
 
 from .orb_cy import _orb_loop
 
@@ -166,6 +167,7 @@ class ORB(FeatureDetector, DescriptorExtractor):
             Input image.
 
         """
+        assert_nD(image, 2)
 
         pyramid = self._build_pyramid(image)
 
@@ -237,6 +239,7 @@ class ORB(FeatureDetector, DescriptorExtractor):
             Corresponding orientations in radians.
 
         """
+        assert_nD(image, 2)
 
         pyramid = self._build_pyramid(image)
 
@@ -282,6 +285,7 @@ class ORB(FeatureDetector, DescriptorExtractor):
             Input image.
 
         """
+        assert_nD(image, 2)
 
         pyramid = self._build_pyramid(image)
 
diff --git a/skimage/feature/template.py b/skimage/feature/template.py
index fbf95866..d7566310 100644
--- a/skimage/feature/template.py
+++ b/skimage/feature/template.py
@@ -2,6 +2,7 @@ import numpy as np
 from scipy.signal import fftconvolve
 
 from skimage.util import pad
+from skimage._shared.utils import assert_nD
 
 
 def _window_sum_2d(image, window_shape):
@@ -102,9 +103,8 @@ def match_template(image, template, pad_input=False, mode='constant',
            [ 0.   ,  0.   ,  0.   ,  0.125, -1.   ,  0.125],
            [ 0.   ,  0.   ,  0.   ,  0.125,  0.125,  0.125]], dtype=float32)
     """
+    assert_nD(image, (2, 3))
 
-    if image.ndim not in (2, 3) or template.ndim not in (2, 3):
-        raise ValueError("Only 2- and 3-D images supported.")
     if image.ndim < template.ndim:
         raise ValueError("Dimensionality of template must be less than or "
                          "equal to the dimensionality of image.")
diff --git a/skimage/feature/tests/test_canny.py b/skimage/feature/tests/test_canny.py
index 43db5037..400a7f70 100644
--- a/skimage/feature/tests/test_canny.py
+++ b/skimage/feature/tests/test_canny.py
@@ -60,7 +60,7 @@ class TestCanny(unittest.TestCase):
         self.assertTrue(point_count < 1600)
 
     def test_image_shape(self):
-        self.assertRaises(TypeError, F.canny, np.zeros((20, 20, 20)), 4, 0, 0)
+        self.assertRaises(ValueError, F.canny, np.zeros((20, 20, 20)), 4, 0, 0)
 
     def test_mask_none(self):
         result1 = F.canny(np.zeros((20, 20)), 4, 0, 0, np.ones((20, 20), bool))
diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py
index fa705f49..b1e1024a 100644
--- a/skimage/feature/texture.py
+++ b/skimage/feature/texture.py
@@ -3,7 +3,7 @@ Methods to characterize image textures.
 """
 
 import numpy as np
-
+from skimage._shared.utils import assert_nD
 from ._texture import _glcm_loop, _local_binary_pattern
 
 
@@ -89,17 +89,17 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False,
            [0, 0, 0, 0]], dtype=uint32)
 
     """
+    assert_nD(image, 2)
+    assert_nD(distances, 1, 'distances')
+    assert_nD(angles, 1, 'angles')
 
     assert levels <= 256
     image = np.ascontiguousarray(image)
-    assert image.ndim == 2
     assert image.min() >= 0
     assert image.max() < levels
     image = image.astype(np.uint8)
     distances = np.ascontiguousarray(distances, dtype=np.float64)
     angles = np.ascontiguousarray(angles, dtype=np.float64)
-    assert distances.ndim == 1
-    assert angles.ndim == 1
 
     P = np.zeros((levels, levels, len(distances), len(angles)),
                  dtype=np.uint32, order='C')
@@ -179,8 +179,8 @@ def greycoprops(P, prop='contrast'):
            [ 1.25      ,  2.75      ]])
 
     """
+    assert_nD(P, 4, 'P')
 
-    assert P.ndim == 4
     (num_level, num_level2, num_dist, num_angle) = P.shape
     assert num_level == num_level2
     assert num_dist > 0
@@ -279,6 +279,7 @@ def local_binary_pattern(image, P, R, method='default'):
            http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.214.6851,
            2004.
     """
+    assert_nD(image, 2)
 
     methods = {
         'default': ord('D'),
diff --git a/skimage/feature/util.py b/skimage/feature/util.py
index 5a3e5687..e70be719 100644
--- a/skimage/feature/util.py
+++ b/skimage/feature/util.py
@@ -1,6 +1,7 @@
 import numpy as np
 
 from skimage.util import img_as_float
+from skimage._shared.utils import assert_nD
 
 
 class FeatureDetector(object):
@@ -124,9 +125,7 @@ def plot_matches(ax, image1, image2, keypoints1, keypoints2, matches,
 
 def _prepare_grayscale_input_2D(image):
     image = np.squeeze(image)
-    if image.ndim != 2:
-        raise ValueError("Only 2-D gray-scale images supported.")
-
+    assert_nD(image, 2)
     return img_as_float(image)
 
 
diff --git a/skimage/filter/_gabor.py b/skimage/filter/_gabor.py
index f766ac74..b1f3fe3b 100644
--- a/skimage/filter/_gabor.py
+++ b/skimage/filter/_gabor.py
@@ -1,5 +1,6 @@
 import numpy as np
 from scipy import ndimage
+from skimage._shared.utils import assert_nD
 
 
 __all__ = ['gabor_kernel', 'gabor_filter']
@@ -112,7 +113,7 @@ def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None,
     .. [2] http://mplab.ucsd.edu/tutorials/gabor.pdf
 
     """
-
+    assert_nD(image, 2)
     g = gabor_kernel(frequency, theta, bandwidth, sigma_x, sigma_y, offset)
 
     filtered_real = ndimage.convolve(image, np.real(g), mode=mode, cval=cval)
diff --git a/skimage/filter/edges.py b/skimage/filter/edges.py
index 945c9ee1..ff3b6b03 100644
--- a/skimage/filter/edges.py
+++ b/skimage/filter/edges.py
@@ -81,7 +81,7 @@ def sobel(image, mask=None):
     Note that ``scipy.ndimage.sobel`` returns a directional Sobel which
     has to be further processed to perform edge detection.
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     return np.sqrt(hsobel(image, mask)**2 + vsobel(image, mask)**2)
 
 
@@ -112,7 +112,7 @@ def hsobel(image, mask=None):
      -1  -2  -1
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, HSOBEL_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -145,7 +145,7 @@ def vsobel(image, mask=None):
       1   0  -1
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, VSOBEL_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -215,7 +215,7 @@ def hscharr(image, mask=None):
            of Kernel Based Image Derivatives.
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, HSCHARR_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -253,7 +253,7 @@ def vscharr(image, mask=None):
            of Kernel Based Image Derivatives.
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, VSCHARR_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -281,7 +281,7 @@ def prewitt(image, mask=None):
     Return the square root of the sum of squares of the horizontal
     and vertical Prewitt transforms.
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     return np.sqrt(hprewitt(image, mask)**2 + vprewitt(image, mask)**2)
 
 
@@ -312,7 +312,7 @@ def hprewitt(image, mask=None):
      -1  -1  -1
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, HPREWITT_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -345,7 +345,7 @@ def vprewitt(image, mask=None):
       1   0  -1
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, VPREWITT_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -368,7 +368,7 @@ def roberts(image, mask=None):
     output : 2-D array
         The Roberts' Cross edge map.
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     return np.sqrt(roberts_positive_diagonal(image, mask)**2 +
                    roberts_negative_diagonal(image, mask)**2)
 
@@ -402,7 +402,7 @@ def roberts_positive_diagonal(image, mask=None):
       0  -1
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, ROBERTS_PD_WEIGHTS))
     return _mask_filter_result(result, mask)
@@ -437,7 +437,7 @@ def roberts_negative_diagonal(image, mask=None):
      -1   0
 
     """
-    assert_nD(image)
+    assert_nD(image, 2)
     image = img_as_float(image)
     result = np.abs(convolve(image, ROBERTS_ND_WEIGHTS))
     return _mask_filter_result(result, mask)
diff --git a/skimage/filter/lpi_filter.py b/skimage/filter/lpi_filter.py
index ef85c5cf..b02d4eb1 100644
--- a/skimage/filter/lpi_filter.py
+++ b/skimage/filter/lpi_filter.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 from scipy.fftpack import ifftshift
+from skimage._shared.utils import assert_nD
 
 eps = np.finfo(float).eps
 
@@ -118,6 +119,7 @@ class LPIFilter2D(object):
         data : (M,N) ndarray
 
         """
+        assert_nD(data, 2, 'data')
         F, G = self._prepare(data)
         out = np.dual.ifftn(F * G)
         out = np.abs(_centre(out, data.shape))
@@ -155,6 +157,7 @@ def forward(data, impulse_response=None, filter_params={},
     >>> filtered = forward(data.coins(), filt_func)
 
     """
+    assert_nD(data, 2, 'data')
     if predefined_filter is None:
         predefined_filter = LPIFilter2D(impulse_response, **filter_params)
     return predefined_filter(data)
@@ -184,6 +187,7 @@ def inverse(data, impulse_response=None, filter_params={}, max_gain=2,
         images, construct the LPIFilter2D and specify it here.
 
     """
+    assert_nD(data, 2, 'data')
     if predefined_filter is None:
         filt = LPIFilter2D(impulse_response, **filter_params)
     else:
@@ -222,6 +226,11 @@ def wiener(data, impulse_response=None, filter_params={}, K=0.25,
         images, construct the LPIFilter2D and specify it here.
 
     """
+    assert_nD(data, 2, 'data')
+
+    if not isinstance(K, float):
+        assert_nD(K, 2, 'K')
+
     if predefined_filter is None:
         filt = LPIFilter2D(impulse_response, **filter_params)
     else:
diff --git a/skimage/filter/rank/_percentile.py b/skimage/filter/rank/_percentile.py
index 9b93195b..e57f5df1 100644
--- a/skimage/filter/rank/_percentile.py
+++ b/skimage/filter/rank/_percentile.py
@@ -23,6 +23,7 @@ References
 """
 
 import numpy as np
+from skimage._shared.utils import assert_nD
 
 from . import percentile_cy
 from .generic import _handle_input
@@ -37,6 +38,7 @@ __all__ = ['autolevel_percentile', 'gradient_percentile',
 def _apply(func, image, selem, out, mask, shift_x, shift_y, p0, p1,
            out_dtype=None):
 
+    assert_nD(image, 2)
     image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask,
                                                      out_dtype)
 
diff --git a/skimage/filter/rank/bilateral.py b/skimage/filter/rank/bilateral.py
index aeb318d1..05a5664a 100644
--- a/skimage/filter/rank/bilateral.py
+++ b/skimage/filter/rank/bilateral.py
@@ -25,6 +25,7 @@ References
 
 import numpy as np
 from skimage import img_as_ubyte
+from skimage._shared.utils import assert_nD
 
 from . import bilateral_cy
 from .generic import _handle_input
@@ -36,6 +37,7 @@ __all__ = ['mean_bilateral', 'pop_bilateral', 'sum_bilateral']
 def _apply(func, image, selem, out, mask, shift_x, shift_y, s0, s1,
            out_dtype=None):
 
+    assert_nD(image, 2)
     image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask,
                                                      out_dtype)
 
diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py
index a474d050..d0a9baec 100644
--- a/skimage/filter/rank/generic.py
+++ b/skimage/filter/rank/generic.py
@@ -19,6 +19,7 @@ References
 import warnings
 import numpy as np
 from skimage import img_as_ubyte
+from skimage._shared.utils import assert_nD
 
 from . import generic_cy
 
@@ -30,6 +31,7 @@ __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean',
 
 def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1):
 
+    assert_nD(image, 2)
     if image.dtype not in (np.uint8, np.uint16):
         image = img_as_ubyte(image)
 
diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py
index 3985e278..5a737293 100644
--- a/skimage/filter/thresholding.py
+++ b/skimage/filter/thresholding.py
@@ -6,6 +6,7 @@ __all__ = ['threshold_adaptive',
 import numpy as np
 import scipy.ndimage
 from skimage.exposure import histogram
+from skimage._shared.utils import assert_nD
 
 
 def threshold_adaptive(image, block_size, method='gaussian', offset=0,
@@ -65,6 +66,7 @@ def threshold_adaptive(image, block_size, method='gaussian', offset=0,
     >>> func = lambda arr: arr.mean()
     >>> binary_image2 = threshold_adaptive(image, 15, 'generic', param=func)
     """
+    assert_nD(image, 2)
     thresh_image = np.zeros(image.shape, 'double')
     if method == 'generic':
         scipy.ndimage.generic_filter(image, param, block_size,