diff --git a/doc/examples/applications/plot_coins_segmentation.py b/doc/examples/applications/plot_coins_segmentation.py
index 98f5212f..fe165760 100644
--- a/doc/examples/applications/plot_coins_segmentation.py
+++ b/doc/examples/applications/plot_coins_segmentation.py
@@ -71,9 +71,9 @@ ax.set_title('Canny detector')
 These contours are then filled using mathematical morphology.
 """
 
-from scipy import ndimage
+from scipy import ndimage as ndi
 
-fill_coins = ndimage.binary_fill_holes(edges)
+fill_coins = ndi.binary_fill_holes(edges)
 
 fig, ax = plt.subplots(figsize=(4, 3))
 ax.imshow(fill_coins, cmap=plt.cm.gray, interpolation='nearest')
@@ -158,8 +158,8 @@ individually.
 
 from skimage.color import label2rgb
 
-segmentation = ndimage.binary_fill_holes(segmentation - 1)
-labeled_coins, _ = ndimage.label(segmentation)
+segmentation = ndi.binary_fill_holes(segmentation - 1)
+labeled_coins, _ = ndi.label(segmentation)
 image_label_overlay = label2rgb(labeled_coins, image=coins)
 
 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3))
diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py
index 4defb8e6..e2dbd71b 100644
--- a/doc/examples/applications/plot_rank_filters.py
+++ b/doc/examples/applications/plot_rank_filters.py
@@ -547,7 +547,7 @@ available in `skimage`.
 
 from time import time
 
-from scipy.ndimage.filters import percentile_filter
+from scipy.ndimage import percentile_filter
 from skimage.morphology import dilation
 from skimage.filters.rank import median, maximum
 
diff --git a/doc/examples/plot_canny.py b/doc/examples/plot_canny.py
index 82a930e5..f01e6e2d 100644
--- a/doc/examples/plot_canny.py
+++ b/doc/examples/plot_canny.py
@@ -17,7 +17,7 @@ the hysteresis thresholding.
 """
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy import ndimage
+from scipy import ndimage as ndi
 
 from skimage import feature
 
@@ -26,8 +26,8 @@ from skimage import feature
 im = np.zeros((128, 128))
 im[32:-32, 32:-32] = 1
 
-im = ndimage.rotate(im, 15, mode='constant')
-im = ndimage.gaussian_filter(im, 4)
+im = ndi.rotate(im, 15, mode='constant')
+im = ndi.gaussian_filter(im, 4)
 im += 0.2 * np.random.random(im.shape)
 
 # Compute the Canny filter for two values of sigma
diff --git a/doc/examples/plot_gabor.py b/doc/examples/plot_gabor.py
index e5cc7760..36d20836 100644
--- a/doc/examples/plot_gabor.py
+++ b/doc/examples/plot_gabor.py
@@ -14,10 +14,9 @@ for classification, which is based on the least squared error for simplicity.
 """
 from __future__ import print_function
 
-import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 
 from skimage import data
 from skimage.util import img_as_float
@@ -27,7 +26,7 @@ from skimage.filters import gabor_kernel
 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')
+        filtered = ndi.convolve(image, kernel, mode='wrap')
         feats[k, 0] = filtered.mean()
         feats[k, 1] = filtered.var()
     return feats
@@ -71,23 +70,23 @@ ref_feats[2, :, :] = compute_feats(wall, kernels)
 print('Rotated images matched against references using Gabor filter banks:')
 
 print('original: brick, rotated: 30deg, match result: ', end='')
-feats = compute_feats(nd.rotate(brick, angle=190, reshape=False), kernels)
+feats = compute_feats(ndi.rotate(brick, angle=190, reshape=False), kernels)
 print(image_names[match(feats, ref_feats)])
 
 print('original: brick, rotated: 70deg, match result: ', end='')
-feats = compute_feats(nd.rotate(brick, angle=70, reshape=False), kernels)
+feats = compute_feats(ndi.rotate(brick, angle=70, reshape=False), kernels)
 print(image_names[match(feats, ref_feats)])
 
 print('original: grass, rotated: 145deg, match result: ', end='')
-feats = compute_feats(nd.rotate(grass, angle=145, reshape=False), kernels)
+feats = compute_feats(ndi.rotate(grass, angle=145, reshape=False), kernels)
 print(image_names[match(feats, ref_feats)])
 
 
 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)
+    return np.sqrt(ndi.convolve(image, np.real(kernel), mode='wrap')**2 +
+                   ndi.convolve(image, np.imag(kernel), mode='wrap')**2)
 
 # Plot a selection of the filter bank kernels and their responses.
 results = []
diff --git a/doc/examples/plot_join_segmentations.py b/doc/examples/plot_join_segmentations.py
index c10105c5..5fded86e 100644
--- a/doc/examples/plot_join_segmentations.py
+++ b/doc/examples/plot_join_segmentations.py
@@ -11,7 +11,7 @@ segmentations.
 
 """
 import numpy as np
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 import matplotlib.pyplot as plt
 
 from skimage.filters import sobel
@@ -30,7 +30,7 @@ markers[coins < 30.0 / 255] = background
 markers[coins > 150.0 / 255] = foreground
 
 ws = watershed(edges, markers)
-seg1 = nd.label(ws == foreground)[0]
+seg1 = ndi.label(ws == foreground)[0]
 
 # make segmentation using SLIC superpixels
 seg2 = slic(coins, n_segments=117, max_iter=160, sigma=1, compactness=0.75,
diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/plot_marked_watershed.py
index 766a1d25..234e3818 100644
--- a/doc/examples/plot_marked_watershed.py
+++ b/doc/examples/plot_marked_watershed.py
@@ -14,7 +14,7 @@ See Wikipedia_ for more details on the algorithm.
 
 """
 
-from scipy import ndimage
+from scipy import ndimage as ndi
 import matplotlib.pyplot as plt
 
 from skimage.morphology import watershed, disk
@@ -30,7 +30,7 @@ denoised = rank.median(image, disk(2))
 
 # find continuous region (low gradient) --> markers
 markers = rank.gradient(denoised, disk(5)) < 10
-markers = ndimage.label(markers)[0]
+markers = ndi.label(markers)[0]
 
 #local gradient
 gradient = rank.gradient(denoised, disk(2))
diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py
index d2314d3d..436e439e 100644
--- a/doc/examples/plot_medial_transform.py
+++ b/doc/examples/plot_medial_transform.py
@@ -21,7 +21,7 @@ a skeleton by iterative morphological thinnings.
 
 """
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from skimage.morphology import medial_axis
 import matplotlib.pyplot as plt
 
@@ -43,7 +43,7 @@ def microstructure(l=256):
     generator = np.random.RandomState(1)
     points = l * generator.rand(2, n**2)
     mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
-    mask = ndimage.gaussian_filter(mask, sigma=l/(4.*n))
+    mask = ndi.gaussian_filter(mask, sigma=l/(4.*n))
     return mask > mask.mean()
 
 data = microstructure(l=64)
diff --git a/doc/examples/plot_peak_local_max.py b/doc/examples/plot_peak_local_max.py
index 03981578..aba2aab0 100644
--- a/doc/examples/plot_peak_local_max.py
+++ b/doc/examples/plot_peak_local_max.py
@@ -10,7 +10,7 @@ size of the dilation. Locations where the original image is equal to the
 dilated image are returned as local maxima.
 
 """
-from scipy import ndimage
+from scipy import ndimage as ndi
 import matplotlib.pyplot as plt
 from skimage.feature import peak_local_max
 from skimage import data, img_as_float
@@ -19,7 +19,7 @@ im = img_as_float(data.coins())
 
 # image_max is the dilation of im with a 20*20 structuring element
 # It is used within peak_local_max function
-image_max = ndimage.maximum_filter(im, size=20, mode='constant')
+image_max = ndi.maximum_filter(im, size=20, mode='constant')
 
 # Comparison between image_max and im to find the coordinates of local maxima
 coordinates = peak_local_max(im, min_distance=20)
diff --git a/doc/examples/plot_random_walker_segmentation.py b/doc/examples/plot_random_walker_segmentation.py
index 0c494de5..446e938f 100644
--- a/doc/examples/plot_random_walker_segmentation.py
+++ b/doc/examples/plot_random_walker_segmentation.py
@@ -21,7 +21,6 @@ values, and use the random walker for the segmentation.
 
 """
 import numpy as np
-from scipy import ndimage
 import matplotlib.pyplot as plt
 
 from skimage.segmentation import random_walker
diff --git a/doc/examples/plot_register_translation.py b/doc/examples/plot_register_translation.py
index ec254e16..05cf78bf 100644
--- a/doc/examples/plot_register_translation.py
+++ b/doc/examples/plot_register_translation.py
@@ -21,7 +21,7 @@ import matplotlib.pyplot as plt
 from skimage import data
 from skimage.feature import register_translation
 from skimage.feature.register_translation import _upsampled_dft
-from scipy.ndimage.fourier import fourier_shift
+from scipy.ndimage import fourier_shift
 
 image = data.camera()
 shift = (-2.4, 1.32)
diff --git a/doc/examples/plot_watershed.py b/doc/examples/plot_watershed.py
index 27f7016e..9ff16354 100644
--- a/doc/examples/plot_watershed.py
+++ b/doc/examples/plot_watershed.py
@@ -26,7 +26,7 @@ See Wikipedia_ for more details on the algorithm.
 """
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy import ndimage
+from scipy import ndimage as ndi
 
 from skimage.morphology import watershed
 from skimage.feature import peak_local_max
@@ -42,10 +42,10 @@ image = np.logical_or(mask_circle1, mask_circle2)
 
 # Now we want to separate the two objects in image
 # Generate the markers as local maxima of the distance to the background
-distance = ndimage.distance_transform_edt(image)
+distance = ndi.distance_transform_edt(image)
 local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)),
                             labels=image)
-markers = ndimage.label(local_maxi)[0]
+markers = ndi.label(local_maxi)[0]
 labels = watershed(-distance, markers, mask=image)
 
 fig, axes = plt.subplots(ncols=3, figsize=(8, 2.7))
diff --git a/doc/source/user_guide/numpy_images.txt b/doc/source/user_guide/numpy_images.txt
index 40084f48..e976ae7e 100644
--- a/doc/source/user_guide/numpy_images.txt
+++ b/doc/source/user_guide/numpy_images.txt
@@ -169,8 +169,8 @@ Many functions in scikit-image operate on 3D images directly:
 
     >>> im3d = np.random.rand(100, 1000, 1000)
     >>> from skimage import morphology
-    >>> from scipy import ndimage as nd
-    >>> seeds = nd.label(im3d < 0.1)[0]
+    >>> from scipy import ndimage as ndi
+    >>> seeds = ndi.label(im3d < 0.1)[0]
     >>> ws = morphology.watershed(im3d, seeds)
 
 In many cases,
diff --git a/doc/source/user_guide/tutorial_segmentation.txt b/doc/source/user_guide/tutorial_segmentation.txt
index 2cf45374..76af842e 100644
--- a/doc/source/user_guide/tutorial_segmentation.txt
+++ b/doc/source/user_guide/tutorial_segmentation.txt
@@ -50,8 +50,8 @@ boundary of the coins, or inside the coins.
 
 ::
 
-    >>> from scipy import ndimage
-    >>> fill_coins = ndimage.binary_fill_holes(edges)
+    >>> from scipy import ndimage as ndi
+    >>> fill_coins = ndi.binary_fill_holes(edges)
 
 .. image:: ../auto_examples/applications/images/plot_coins_segmentation_3.png
    :target: ../auto_examples/applications/plot_coins_segmentation.html
@@ -59,7 +59,7 @@ boundary of the coins, or inside the coins.
 
 Now that we have contours that delineate the outer boundary of the coins,
 we fill the inner part of the coins using the
-``ndimage.binary_fill_holes`` function, which uses mathematical morphology
+``ndi.binary_fill_holes`` function, which uses mathematical morphology
 to fill the holes.
 
 .. image:: ../auto_examples/applications/images/plot_coins_segmentation_4.png
@@ -67,12 +67,12 @@ to fill the holes.
    :align: center
 
 Most coins are well segmented out of the background. Small objects from
-the background can be easily removed using the ``ndimage.label``
+the background can be easily removed using the ``ndi.label``
 function to remove objects smaller than a small threshold.
 
 ::
 
-    >>> label_objects, nb_labels = ndimage.label(fill_coins)
+    >>> label_objects, nb_labels = ndi.label(fill_coins)
     >>> sizes = np.bincount(label_objects.ravel())
     >>> mask_sizes = sizes > 20
     >>> mask_sizes[0] = 0
@@ -159,11 +159,11 @@ background.
 
 We remove a few small holes with mathematical morphology::
 
-    >>> segmentation = ndimage.binary_fill_holes(segmentation - 1)
+    >>> segmentation = ndi.binary_fill_holes(segmentation - 1)
 
-We can now label all the coins one by one using ``ndimage.label``::
+We can now label all the coins one by one using ``ndi.label``::
 
-    >>> labeled_coins, _ = ndimage.label(segmentation)
+    >>> labeled_coins, _ = ndi.label(segmentation)
 
 .. image:: ../auto_examples/applications/images/plot_coins_segmentation_9.png
    :target: ../auto_examples/applications/plot_coins_segmentation.html
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 1d4e8339..88532c6a 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -1,6 +1,6 @@
 
 import numpy as np
-from scipy.ndimage.filters import gaussian_filter, gaussian_laplace
+from scipy.ndimage import gaussian_filter, gaussian_laplace
 import itertools as itt
 import math
 from math import sqrt, hypot, log
diff --git a/skimage/feature/brief.py b/skimage/feature/brief.py
index 3a53f89f..ce96afe6 100644
--- a/skimage/feature/brief.py
+++ b/skimage/feature/brief.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy.ndimage.filters import gaussian_filter
+from scipy.ndimage import gaussian_filter
 
 from .util import (DescriptorExtractor, _mask_border_keypoints,
                    _prepare_grayscale_input_2D)
diff --git a/skimage/feature/censure.py b/skimage/feature/censure.py
index d75c03ca..75a7dfb9 100644
--- a/skimage/feature/censure.py
+++ b/skimage/feature/censure.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy.ndimage.filters import maximum_filter, minimum_filter, convolve
+from scipy.ndimage import maximum_filter, minimum_filter, convolve
 
 from ..transform import integral_image
 from ..feature import structure_tensor
diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index 5e2251dc..a8fa7d78 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from scipy import stats
 
 from ..util import img_as_float, pad
@@ -33,8 +33,8 @@ def _compute_derivatives(image, mode='constant', cval=0):
 
     """
 
-    imy = ndimage.sobel(image, axis=0, mode=mode, cval=cval)
-    imx = ndimage.sobel(image, axis=1, mode=mode, cval=cval)
+    imy = ndi.sobel(image, axis=0, mode=mode, cval=cval)
+    imx = ndi.sobel(image, axis=1, mode=mode, cval=cval)
 
     return imx, imy
 
@@ -92,9 +92,9 @@ def structure_tensor(image, sigma=1, mode='constant', cval=0):
     imx, imy = _compute_derivatives(image, mode=mode, cval=cval)
 
     # structure tensore
-    Axx = ndimage.gaussian_filter(imx * imx, sigma, mode=mode, cval=cval)
-    Axy = ndimage.gaussian_filter(imx * imy, sigma, mode=mode, cval=cval)
-    Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode=mode, cval=cval)
+    Axx = ndi.gaussian_filter(imx * imx, sigma, mode=mode, cval=cval)
+    Axy = ndi.gaussian_filter(imx * imy, sigma, mode=mode, cval=cval)
+    Ayy = ndi.gaussian_filter(imy * imy, sigma, mode=mode, cval=cval)
 
     return Axx, Axy, Ayy
 
@@ -165,9 +165,9 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0):
     kernel_xy /= kernel_xx.sum()
     kernel_yy = kernel_xx.transpose()
 
-    Hxx = ndimage.convolve(image, kernel_xx, mode=mode, cval=cval)
-    Hxy = ndimage.convolve(image, kernel_xy, mode=mode, cval=cval)
-    Hyy = ndimage.convolve(image, kernel_yy, mode=mode, cval=cval)
+    Hxx = ndi.convolve(image, kernel_xx, mode=mode, cval=cval)
+    Hxy = ndi.convolve(image, kernel_xy, mode=mode, cval=cval)
+    Hyy = ndi.convolve(image, kernel_yy, mode=mode, cval=cval)
 
     return Hxx, Hxy, Hyy
 
diff --git a/skimage/feature/tests/test_hog.py b/skimage/feature/tests/test_hog.py
index 855bb557..6112d726 100644
--- a/skimage/feature/tests/test_hog.py
+++ b/skimage/feature/tests/test_hog.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from skimage import data
 from skimage import feature
 from skimage import img_as_float
@@ -137,7 +137,7 @@ def test_hog_orientations_circle():
     image = np.zeros((height, width))
     rr, cc = draw.circle(int(height / 2), int(width / 2), int(width / 3))
     image[rr, cc] = 100
-    image = ndimage.gaussian_filter(image, 2)
+    image = ndi.gaussian_filter(image, 2)
 
     for orientations in range(2, 15):
         (hog, hog_img) = feature.hog(image, orientations=orientations,
diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py
index 09e52495..90fe9b62 100644
--- a/skimage/feature/tests/test_peak.py
+++ b/skimage/feature/tests/test_peak.py
@@ -1,7 +1,7 @@
 import numpy as np
 from numpy.testing import (assert_array_almost_equal as assert_close,
                            assert_equal)
-import scipy.ndimage
+from scipy import ndimage as ndi
 from skimage.feature import peak
 
 
@@ -101,8 +101,8 @@ def test_reorder_labels():
     expected = np.zeros(image.shape, float)
     for imin, imax in ((0, 20), (20, 40)):
         for jmin, jmax in ((0, 30), (30, 60)):
-            expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter(
-                image[imin:imax, jmin:jmax], footprint=footprint)
+            expected[imin:imax, jmin:jmax] = ndi.maximum_filter(
+                            image[imin:imax, jmin:jmax], footprint=footprint)
     expected = (expected == image)
     result = peak.peak_local_max(image, labels=labels, min_distance=1,
                                  threshold_rel=0, footprint=footprint,
@@ -119,7 +119,7 @@ def test_indices_with_labels():
     expected = np.zeros(image.shape, float)
     for imin, imax in ((0, 20), (20, 40)):
         for jmin, jmax in ((0, 30), (30, 60)):
-            expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter(
+            expected[imin:imax, jmin:jmax] = ndi.maximum_filter(
                 image[imin:imax, jmin:jmax], footprint=footprint)
     expected = (expected == image)
     result = peak.peak_local_max(image, labels=labels, min_distance=1,
@@ -251,7 +251,7 @@ def test_four_quadrants():
     expected = np.zeros(image.shape, float)
     for imin, imax in ((0, 20), (20, 40)):
         for jmin, jmax in ((0, 30), (30, 60)):
-            expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter(
+            expected[imin:imax, jmin:jmax] = ndi.maximum_filter(
                 image[imin:imax, jmin:jmax], footprint=footprint)
     expected = (expected == image)
     result = peak.peak_local_max(image, labels=labels, footprint=footprint,
diff --git a/skimage/feature/tests/test_register_translation.py b/skimage/feature/tests/test_register_translation.py
index 6494630e..b3a70ad1 100644
--- a/skimage/feature/tests/test_register_translation.py
+++ b/skimage/feature/tests/test_register_translation.py
@@ -4,7 +4,7 @@ from numpy.testing import assert_allclose, assert_raises
 from skimage.feature.register_translation import (register_translation,
                                                   _upsampled_dft)
 from skimage.data import camera
-from scipy.ndimage.fourier import fourier_shift
+from scipy.ndimage import fourier_shift
 
 
 def test_correlation():
diff --git a/skimage/filters/_gabor.py b/skimage/filters/_gabor.py
index 862a593c..3eeccc0f 100644
--- a/skimage/filters/_gabor.py
+++ b/skimage/filters/_gabor.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from .._shared.utils import assert_nD
 
 
@@ -130,10 +130,10 @@ def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None,
     offset : float, optional
         Phase offset of harmonic function in radians.
     mode : string, optional
-        Mode used to convolve image with a kernel, passed to `ndimage.convolve`
+        Mode used to convolve image with a kernel, passed to `ndi.convolve`
     cval : scalar, optional
         Value to fill past edges of input if `mode` of convolution is
-        'constant'. The parameter is passed to `ndimage.convolve`.
+        'constant'. The parameter is passed to `ndi.convolve`.
 
     Returns
     -------
@@ -169,7 +169,7 @@ def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None,
     g = gabor_kernel(frequency, theta, bandwidth, sigma_x, sigma_y, n_stds,
                      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)
+    filtered_real = ndi.convolve(image, np.real(g), mode=mode, cval=cval)
+    filtered_imag = ndi.convolve(image, np.imag(g), mode=mode, cval=cval)
 
     return filtered_real, filtered_imag
diff --git a/skimage/filters/_gaussian.py b/skimage/filters/_gaussian.py
index eedafbdb..6bc12574 100644
--- a/skimage/filters/_gaussian.py
+++ b/skimage/filters/_gaussian.py
@@ -1,6 +1,6 @@
 import collections as coll
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 import warnings
 
 from ..util import img_as_float
@@ -46,7 +46,7 @@ def gaussian_filter(image, sigma, output=None, mode='nearest', cval=0,
 
     Notes
     -----
-    This function is a wrapper around :func:`scipy.ndimage.gaussian_filter`.
+    This function is a wrapper around :func:`scipy.ndi.gaussian_filter`.
 
     Integer arrays are converted to float.
 
@@ -100,4 +100,4 @@ def gaussian_filter(image, sigma, output=None, mode='nearest', cval=0,
         if len(sigma) != image.ndim:
             sigma = np.concatenate((np.asarray(sigma), [0]))
     image = img_as_float(image)
-    return ndimage.gaussian_filter(image, sigma, mode=mode, cval=cval)
+    return ndi.gaussian_filter(image, sigma, mode=mode, cval=cval)
diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py
index 473db40a..d962c562 100644
--- a/skimage/filters/thresholding.py
+++ b/skimage/filters/thresholding.py
@@ -5,7 +5,7 @@ __all__ = ['threshold_adaptive',
            'threshold_li', ]
 
 import numpy as np
-import scipy.ndimage
+from scipy import ndimage as ndi
 from ..exposure import histogram
 from .._shared.utils import assert_nD
 
@@ -70,26 +70,23 @@ def threshold_adaptive(image, block_size, method='gaussian', offset=0,
     assert_nD(image, 2)
     thresh_image = np.zeros(image.shape, 'double')
     if method == 'generic':
-        scipy.ndimage.generic_filter(image, param, block_size,
-                                     output=thresh_image, mode=mode)
+        ndi.generic_filter(image, param, block_size,
+                           output=thresh_image, mode=mode)
     elif method == 'gaussian':
         if param is None:
             # automatically determine sigma which covers > 99% of distribution
             sigma = (block_size - 1) / 6.0
         else:
             sigma = param
-        scipy.ndimage.gaussian_filter(image, sigma, output=thresh_image,
-                                      mode=mode)
+        ndi.gaussian_filter(image, sigma, output=thresh_image, mode=mode)
     elif method == 'mean':
         mask = 1. / block_size * np.ones((block_size,))
         # separation of filters to speedup convolution
-        scipy.ndimage.convolve1d(image, mask, axis=0, output=thresh_image,
-                                 mode=mode)
-        scipy.ndimage.convolve1d(thresh_image, mask, axis=1,
-                                 output=thresh_image, mode=mode)
+        ndi.convolve1d(image, mask, axis=0, output=thresh_image, mode=mode)
+        ndi.convolve1d(thresh_image, mask, axis=1,
+                       output=thresh_image, mode=mode)
     elif method == 'median':
-        scipy.ndimage.median_filter(image, block_size, output=thresh_image,
-                                    mode=mode)
+        ndi.median_filter(image, block_size, output=thresh_image, mode=mode)
 
     return image > (thresh_image - offset)
 
diff --git a/skimage/future/graph/rag.py b/skimage/future/graph/rag.py
index 11441800..2975a027 100644
--- a/skimage/future/graph/rag.py
+++ b/skimage/future/graph/rag.py
@@ -11,8 +11,7 @@ except ImportError:
     warnings.warn(msg)
 
 import numpy as np
-from scipy.ndimage import filters
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 import math
 from ... import draw, measure, segmentation, util, color
 try:
@@ -269,7 +268,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance',
     # The footprint is constructed in such a way that the first
     # element in the array being passed to _add_edge_filter is
     # the central value.
-    fp = nd.generate_binary_structure(labels.ndim, connectivity)
+    fp = ndi.generate_binary_structure(labels.ndim, connectivity)
     for d in range(fp.ndim):
         fp = fp.swapaxes(0, d)
         fp[0, ...] = 0
@@ -286,7 +285,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance',
     #       [0,1,1],
     #       [0,1,1]]
 
-    filters.generic_filter(
+    ndi.generic_filter(
         labels,
         function=_add_edge_filter,
         footprint=fp,
diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py
index ca778272..5f91de3e 100644
--- a/skimage/measure/_regionprops.py
+++ b/skimage/measure/_regionprops.py
@@ -1,8 +1,7 @@
 # coding: utf-8
-import warnings
 from math import sqrt, atan2, pi as PI
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 
 from ._label import label
 from . import _moments
@@ -172,7 +171,7 @@ class _RegionProperties(object):
 
     @_cached_property
     def filled_image(self):
-        return ndimage.binary_fill_holes(self.image, STREL_8)
+        return ndi.binary_fill_holes(self.image, STREL_8)
 
     @_cached_property
     def image(self):
@@ -494,7 +493,7 @@ def regionprops(label_image, intensity_image=None, cache=True):
 
     regions = []
 
-    objects = ndimage.find_objects(label_image)
+    objects = ndi.find_objects(label_image)
     for i, sl in enumerate(objects):
         if sl is None:
             continue
@@ -534,7 +533,7 @@ def perimeter(image, neighbourhood=4):
     else:
         strel = STREL_8
     image = image.astype(np.uint8)
-    eroded_image = ndimage.binary_erosion(image, strel, border_value=0)
+    eroded_image = ndi.binary_erosion(image, strel, border_value=0)
     border_image = image - eroded_image
 
     perimeter_weights = np.zeros(50, dtype=np.double)
@@ -543,10 +542,10 @@ def perimeter(image, neighbourhood=4):
     perimeter_weights[[13, 23]] = (1 + sqrt(2)) / 2
 
 
-    perimeter_image = ndimage.convolve(border_image, np.array([[10, 2, 10],
-                                                               [ 2, 1,  2],
-                                                               [10, 2, 10]]),
-                                       mode='constant', cval=0)
+    perimeter_image = ndi.convolve(border_image, np.array([[10, 2, 10],
+                                                           [ 2, 1,  2],
+                                                           [10, 2, 10]]),
+                                   mode='constant', cval=0)
 
     # You can also write
     # return perimeter_weights[perimeter_image].sum()
diff --git a/skimage/measure/_structural_similarity.py b/skimage/measure/_structural_similarity.py
index bfe02e6a..0b431a21 100644
--- a/skimage/measure/_structural_similarity.py
+++ b/skimage/measure/_structural_similarity.py
@@ -3,7 +3,7 @@ from __future__ import division
 __all__ = ['structural_similarity']
 
 import numpy as np
-from scipy.ndimage.filters import uniform_filter, gaussian_filter
+from scipy.ndimage import uniform_filter, gaussian_filter
 
 from ..util.dtype import dtype_range
 from ..util.arraypad import crop
diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py
index 15353e51..819b8cff 100644
--- a/skimage/measure/profile.py
+++ b/skimage/measure/profile.py
@@ -1,5 +1,5 @@
 import numpy as np
-import scipy.ndimage as nd
+from scipy import ndimage as ndi
 
 
 def profile_line(img, src, dst, linewidth=1,
@@ -52,13 +52,13 @@ def profile_line(img, src, dst, linewidth=1,
     """
     perp_lines = _line_profile_coordinates(src, dst, linewidth=linewidth)
     if img.ndim == 3:
-        pixels = [nd.map_coordinates(img[..., i], perp_lines,
-                                     order=order, mode=mode, cval=cval)
+        pixels = [ndi.map_coordinates(img[..., i], perp_lines,
+                                      order=order, mode=mode, cval=cval)
                   for i in range(img.shape[2])]
         pixels = np.transpose(np.asarray(pixels), (1, 2, 0))
     else:
-        pixels = nd.map_coordinates(img, perp_lines,
-                                    order=order, mode=mode, cval=cval)
+        pixels = ndi.map_coordinates(img, perp_lines,
+                                     order=order, mode=mode, cval=cval)
     intensities = pixels.mean(axis=1)
 
     return intensities
diff --git a/skimage/morphology/_skeletonize.py b/skimage/morphology/_skeletonize.py
index d310bad3..e892c185 100644
--- a/skimage/morphology/_skeletonize.py
+++ b/skimage/morphology/_skeletonize.py
@@ -3,7 +3,7 @@ Algorithms for computing the skeleton of a binary image
 """
 
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 
 from ._skeletonize_cy import _skeletonize_loop, _table_lookup_index
 
@@ -119,7 +119,7 @@ def skeletonize(image):
         pixel_removed = False
 
         # assign each pixel a unique value based on its foreground neighbours
-        neighbours = ndimage.correlate(skeleton, mask, mode='constant')
+        neighbours = ndi.correlate(skeleton, mask, mode='constant')
 
         # ignore background
         neighbours *= skeleton
@@ -138,7 +138,7 @@ def skeletonize(image):
             skeleton[code_mask] = 0
 
         # pass 2 - remove the 2's and 3's
-        neighbours = ndimage.correlate(skeleton, mask, mode='constant')
+        neighbours = ndi.correlate(skeleton, mask, mode='constant')
         neighbours *= skeleton
         codes = np.take(lut, neighbours)
         code_mask = (codes == 2)
@@ -154,7 +154,7 @@ def skeletonize(image):
 
 # --------- Skeletonization by medial axis transform --------
 
-_eight_connect = ndimage.generate_binary_structure(2, 2)
+_eight_connect = ndi.generate_binary_structure(2, 2)
 
 
 def medial_axis(image, mask=None, return_distance=False):
@@ -247,17 +247,17 @@ def medial_axis(image, mask=None, return_distance=False):
     center_is_foreground = (np.arange(512) & 2**4).astype(bool)
     table = (center_is_foreground  # condition 1.
                 &
-            (np.array([ndimage.label(_pattern_of(index), _eight_connect)[1] !=
-                        ndimage.label(_pattern_of(index & ~ 2**4),
+            (np.array([ndi.label(_pattern_of(index), _eight_connect)[1] !=
+                       ndi.label(_pattern_of(index & ~ 2**4),
                                     _eight_connect)[1]
-                        for index in range(512)])  # condition 2
+                       for index in range(512)])  # condition 2
                 |
         np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)]))
         # condition 3
             )
 
     # Build distance transform
-    distance = ndimage.distance_transform_edt(masked_image)
+    distance = ndi.distance_transform_edt(masked_image)
     if return_distance:
         store_distance = distance.copy()
 
diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py
index 0ff5c8f7..fffd8993 100644
--- a/skimage/morphology/binary.py
+++ b/skimage/morphology/binary.py
@@ -2,7 +2,7 @@
 Binary morphological operations
 """
 import numpy as np
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 from .misc import default_selem
 
 
@@ -39,7 +39,7 @@ def binary_erosion(image, selem=None, out=None):
     """
     if out is None:
         out = np.empty(image.shape, dtype=np.bool)
-    nd.binary_erosion(image, structure=selem, output=out)
+    ndi.binary_erosion(image, structure=selem, output=out)
     return out
 
 
@@ -74,7 +74,7 @@ def binary_dilation(image, selem=None, out=None):
     """
     if out is None:
         out = np.empty(image.shape, dtype=np.bool)
-    nd.binary_dilation(image, structure=selem, output=out)
+    ndi.binary_dilation(image, structure=selem, output=out)
     return out
 
 
diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py
index fdc8381c..d8d99edf 100644
--- a/skimage/morphology/grey.py
+++ b/skimage/morphology/grey.py
@@ -3,7 +3,7 @@ Grayscale morphological operations
 """
 import functools
 import numpy as np
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 from .misc import default_selem
 from ..util import pad, crop
 
@@ -53,7 +53,7 @@ def _invert_selem(selem):
     """Change the order of the values in `selem`.
 
     This is a patch for the *weird* footprint inversion in
-    `nd.grey_morphology` [1]_.
+    `ndi.grey_morphology` [1]_.
 
     Parameters
     ----------
@@ -182,7 +182,7 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False):
     selem = _shift_selem(selem, shift_x, shift_y)
     if out is None:
         out = np.empty_like(image)
-    nd.grey_erosion(image, footprint=selem, output=out)
+    ndi.grey_erosion(image, footprint=selem, output=out)
     return out
 
 
@@ -243,12 +243,12 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False):
     # Inside ndimage.grey_dilation, the structuring element is inverted,
     # eg. `selem = selem[::-1, ::-1]` for 2D [1]_, for reasons unknown to
     # this author (@jni). To "patch" this behaviour, we invert our own
-    # selem before passing it to `nd.grey_dilation`.
+    # selem before passing it to `ndi.grey_dilation`.
     # [1] https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285
     selem = _invert_selem(selem)
     if out is None:
         out = np.empty_like(image)
-    nd.grey_dilation(image, footprint=selem, output=out)
+    ndi.grey_dilation(image, footprint=selem, output=out)
     return out
 
 
@@ -401,7 +401,7 @@ def white_tophat(image, selem=None, out=None):
         return out
     elif out is None:
         out = np.empty_like(image)
-    out = nd.white_tophat(image, footprint=selem, output=out)
+    out = ndi.white_tophat(image, footprint=selem, output=out)
     return out
 
 
diff --git a/skimage/morphology/misc.py b/skimage/morphology/misc.py
index 663e1e02..e743a398 100644
--- a/skimage/morphology/misc.py
+++ b/skimage/morphology/misc.py
@@ -1,7 +1,7 @@
 import numpy as np
 import functools
 import warnings
-import scipy.ndimage as nd
+from scipy import ndimage as ndi
 from .selem import _default_selem
 
 # Our function names don't exactly correspond to ndimages.
@@ -102,9 +102,9 @@ def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False):
         return out
 
     if out.dtype == bool:
-        selem = nd.generate_binary_structure(ar.ndim, connectivity)
+        selem = ndi.generate_binary_structure(ar.ndim, connectivity)
         ccs = np.zeros_like(ar, dtype=np.int32)
-        nd.label(ar, selem, output=ccs)
+        ndi.label(ar, selem, output=ccs)
     else:
         ccs = out
 
diff --git a/skimage/morphology/selem.py b/skimage/morphology/selem.py
index 937ff45d..53972cdd 100644
--- a/skimage/morphology/selem.py
+++ b/skimage/morphology/selem.py
@@ -4,7 +4,7 @@
 """
 
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from .. import draw
 
 def square(width, dtype=np.uint8):
@@ -360,4 +360,4 @@ def _default_selem(ndim):
         are 1 and 0 otherwise.
 
     """
-    return ndimage.morphology.generate_binary_structure(ndim, 1)
+    return ndi.morphology.generate_binary_structure(ndim, 1)
diff --git a/skimage/morphology/tests/test_binary.py b/skimage/morphology/tests/test_binary.py
index 6d8097b4..b361af64 100644
--- a/skimage/morphology/tests/test_binary.py
+++ b/skimage/morphology/tests/test_binary.py
@@ -4,7 +4,7 @@ from numpy import testing
 from skimage import data, color
 from skimage.util import img_as_bool
 from skimage.morphology import binary, grey, selem
-from scipy import ndimage
+from scipy import ndimage as ndi
 
 
 img = color.rgb2gray(data.astronaut())
@@ -96,7 +96,7 @@ def test_3d_fallback_default_selem():
 
     # expect a "hyper-cross" centered in the 5x5x5:
     image_expected = np.zeros((7, 7, 7), dtype=bool)
-    image_expected[2:5, 2:5, 2:5] = ndimage.generate_binary_structure(3, 1)
+    image_expected[2:5, 2:5, 2:5] = ndi.generate_binary_structure(3, 1)
     testing.assert_array_equal(opened, image_expected)
 
 def test_3d_fallback_cube_selem():
@@ -119,9 +119,9 @@ def test_2d_ndimage_equivalence():
     bin_opened = binary.binary_opening(image)
     bin_closed = binary.binary_closing(image)
 
-    selem = ndimage.generate_binary_structure(2, 1)
-    ndimage_opened = ndimage.binary_opening(image, structure=selem)
-    ndimage_closed = ndimage.binary_closing(image, structure=selem)
+    selem = ndi.generate_binary_structure(2, 1)
+    ndimage_opened = ndi.binary_opening(image, structure=selem)
+    ndimage_closed = ndi.binary_closing(image, structure=selem)
 
     testing.assert_array_equal(bin_opened, ndimage_opened)
     testing.assert_array_equal(bin_closed, ndimage_closed)
diff --git a/skimage/morphology/tests/test_grey.py b/skimage/morphology/tests/test_grey.py
index 576404f2..cbf7d2b8 100644
--- a/skimage/morphology/tests/test_grey.py
+++ b/skimage/morphology/tests/test_grey.py
@@ -2,9 +2,8 @@ import os.path
 
 import numpy as np
 from numpy import testing
-from scipy import ndimage
+from scipy import ndimage as ndi
 
-import skimage
 from skimage import data_dir, img_as_uint
 from skimage.morphology import grey, selem
 from skimage._shared._warnings import expected_warnings
@@ -152,7 +151,7 @@ def test_3d_fallback_default_selem():
 
     # expect a "hyper-cross" centered in the 5x5x5:
     image_expected = np.zeros((7, 7, 7), dtype=bool)
-    image_expected[2:5, 2:5, 2:5] = ndimage.generate_binary_structure(3, 1)
+    image_expected[2:5, 2:5, 2:5] = ndi.generate_binary_structure(3, 1)
     testing.assert_array_equal(opened, image_expected)
 
 def test_3d_fallback_cube_selem():
@@ -174,9 +173,9 @@ def test_3d_fallback_white_tophat():
 
     with expected_warnings(['operator.*deprecated|\A\Z']):
         new_image = grey.white_tophat(image)
-    footprint = ndimage.generate_binary_structure(3,1)
+    footprint = ndi.generate_binary_structure(3,1)
     with expected_warnings(['operator.*deprecated|\A\Z']):
-        image_expected = ndimage.white_tophat(image,footprint=footprint)
+        image_expected = ndi.white_tophat(image,footprint=footprint)
     testing.assert_array_equal(new_image, image_expected)
 
 def test_3d_fallback_black_tophat():
@@ -187,9 +186,9 @@ def test_3d_fallback_black_tophat():
 
     with expected_warnings(['operator.*deprecated|\A\Z']):
         new_image = grey.black_tophat(image)
-    footprint = ndimage.generate_binary_structure(3,1)
+    footprint = ndi.generate_binary_structure(3,1)
     with expected_warnings(['operator.*deprecated|\A\Z']):
-        image_expected = ndimage.black_tophat(image,footprint=footprint)
+        image_expected = ndi.black_tophat(image,footprint=footprint)
     testing.assert_array_equal(new_image, image_expected)
 
 def test_2d_ndimage_equivalence():
@@ -201,9 +200,9 @@ def test_2d_ndimage_equivalence():
     opened = grey.opening(image)
     closed = grey.closing(image)
 
-    selem = ndimage.generate_binary_structure(2, 1)
-    ndimage_opened = ndimage.grey_opening(image, footprint=selem)
-    ndimage_closed = ndimage.grey_closing(image, footprint=selem)
+    selem = ndi.generate_binary_structure(2, 1)
+    ndimage_opened = ndi.grey_opening(image, footprint=selem)
+    ndimage_closed = ndi.grey_closing(image, footprint=selem)
 
     testing.assert_array_equal(opened, ndimage_opened)
     testing.assert_array_equal(closed, ndimage_closed)
diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py
index 39b464cd..cdaed77d 100644
--- a/skimage/morphology/tests/test_watershed.py
+++ b/skimage/morphology/tests/test_watershed.py
@@ -46,7 +46,7 @@ import math
 import unittest
 
 import numpy as np
-import scipy.ndimage
+from scipy import ndimage as ndi
 
 from skimage.morphology.watershed import watershed, _slow_watershed
 
@@ -380,10 +380,9 @@ class TestWatershed(unittest.TestCase):
             markers[x, y] = idx
             idx += 1
 
-        image = scipy.ndimage.gaussian_filter(image, 4)
+        image = ndi.gaussian_filter(image, 4)
         watershed(image, markers, self.eight)
-        scipy.ndimage.watershed_ift(image.astype(np.uint16), markers,
-                                    self.eight)
+        ndi.watershed_ift(image.astype(np.uint16), markers, self.eight)
 
     def test_watershed10(self):
         "watershed 10"
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index d788d0e2..8a39e19e 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -26,9 +26,8 @@ Original author: Lee Kamentsky
 
 from _heapq import heappush, heappop
 import numpy as np
-import scipy.ndimage
+from scipy import ndimage as ndi
 from ..filters import rank_order
-from .._shared.utils import deprecated
 
 from . import _watershed
 
@@ -114,13 +113,13 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
     >>> # Now we want to separate the two objects in image
     >>> # Generate the markers as local maxima of the distance
     >>> # to the background
-    >>> from scipy import ndimage
-    >>> distance = ndimage.distance_transform_edt(image)
+    >>> from scipy import ndimage as ndi
+    >>> distance = ndi.distance_transform_edt(image)
     >>> from skimage.feature import peak_local_max
     >>> local_maxi = peak_local_max(distance, labels=image,
     ...                             footprint=np.ones((3, 3)),
     ...                             indices=False)
-    >>> markers = ndimage.label(local_maxi)[0]
+    >>> markers = ndi.label(local_maxi)[0]
     >>> labels = watershed(-distance, markers, mask=image)
 
     The algorithm works also for 3-D images, and can be used for example to
@@ -128,7 +127,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
     """
 
     if connectivity is None:
-        c_connectivity = scipy.ndimage.generate_binary_structure(image.ndim, 1)
+        c_connectivity = ndi.generate_binary_structure(image.ndim, 1)
     else:
         c_connectivity = np.array(connectivity, bool)
         if c_connectivity.ndim != image.ndim:
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 04a83608..e9d6287c 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -24,10 +24,9 @@ def test_denoise_tv_chambolle_2d():
     denoised_astro = restoration.denoise_tv_chambolle(img, weight=60.0)
     # which dtype?
     assert denoised_astro.dtype in [np.float, np.float32, np.float64]
-    from scipy import ndimage
-    grad = ndimage.morphological_gradient(img, size=((3, 3)))
-    grad_denoised = ndimage.morphological_gradient(
-        denoised_astro, size=((3, 3)))
+    from scipy import ndimage as ndi
+    grad = ndi.morphological_gradient(img, size=((3, 3)))
+    grad_denoised = ndi.morphological_gradient(denoised_astro, size=((3, 3)))
     # test if the total variation has decreased
     assert grad_denoised.dtype == np.float
     assert (np.sqrt((grad_denoised**2).sum())
diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index f248f613..2e322859 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -2,7 +2,7 @@ from os.path import abspath, dirname, join as pjoin
 
 import numpy as np
 from scipy.signal import convolve2d
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 
 import skimage
 from skimage.data import camera
@@ -61,10 +61,10 @@ def test_image_shape():
     """
     point = np.zeros((5, 5), np.float)
     point[2, 2] = 1.
-    psf = nd.gaussian_filter(point, sigma=1.)
+    psf = ndi.gaussian_filter(point, sigma=1.)
     # image shape: (45, 45), as reported in #1172
     image = skimage.img_as_float(camera()[110:155, 225:270]) # just the face
-    image_conv = nd.convolve(image, psf)
+    image_conv = ndi.convolve(image, psf)
     deconv_sup = restoration.wiener(image_conv, psf, 1)
     deconv_un = restoration.unsupervised_wiener(image_conv, psf)[0]
     # test the shape
diff --git a/skimage/segmentation/_felzenszwalb_cy.pyx b/skimage/segmentation/_felzenszwalb_cy.pyx
index 0ab28e48..9e68277b 100644
--- a/skimage/segmentation/_felzenszwalb_cy.pyx
+++ b/skimage/segmentation/_felzenszwalb_cy.pyx
@@ -3,7 +3,7 @@
 #cython: nonecheck=False
 #cython: wraparound=False
 import numpy as np
-import scipy
+from scipy import ndimage as ndi
 
 cimport cython
 cimport numpy as cnp
@@ -47,7 +47,7 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8,
 
     # rescale scale to behave like in reference implementation
     scale = float(scale) / 255.
-    image = scipy.ndimage.gaussian_filter(image, sigma=sigma)
+    image = ndi.gaussian_filter(image, sigma=sigma)
 
     # compute edge weights in 8 connectivity:
     right_cost = np.abs((image[1:, :] - image[:-1, :]))
diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx
index fef74a9e..c3871806 100644
--- a/skimage/segmentation/_quickshift.pyx
+++ b/skimage/segmentation/_quickshift.pyx
@@ -3,7 +3,7 @@
 #cython: nonecheck=False
 #cython: wraparound=False
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from itertools import product
 
 cimport numpy as cnp
@@ -69,7 +69,7 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10,
             ValueError("Only RGB images can be converted to Lab space.")
         image = rgb2lab(image)
 
-    image = ndimage.gaussian_filter(img_as_float(image), [sigma, sigma, 0])
+    image = ndi.gaussian_filter(img_as_float(image), [sigma, sigma, 0])
     cdef cnp.ndarray[dtype=cnp.float_t, ndim=3, mode="c"] image_c \
             = np.ascontiguousarray(image) * ratio
 
diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py
index 9400e306..1a7e1245 100644
--- a/skimage/segmentation/boundaries.py
+++ b/skimage/segmentation/boundaries.py
@@ -1,7 +1,7 @@
 from __future__ import division
 
 import numpy as np
-from scipy import ndimage as nd
+from scipy import ndimage as ndi
 from ..morphology import dilation, erosion, square
 from ..util import img_as_float, view_as_windows, pad
 from ..color import gray2rgb
@@ -146,7 +146,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0):
            [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
     """
     ndim = label_img.ndim
-    selem = nd.generate_binary_structure(ndim, connectivity)
+    selem = ndi.generate_binary_structure(ndim, connectivity)
     if mode != 'subpixel':
         boundaries = dilation(label_img, selem) != erosion(label_img, selem)
         if mode == 'inner':
@@ -155,7 +155,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0):
         elif mode == 'outer':
             max_label = np.iinfo(label_img.dtype).max
             background_image = (label_img == background)
-            selem = nd.generate_binary_structure(ndim, ndim)
+            selem = ndi.generate_binary_structure(ndim, ndim)
             inverted_background = np.array(label_img, copy=True)
             inverted_background[background_image] = max_label
             adjacent_objects = ((dilation(label_img, selem) !=
@@ -205,10 +205,10 @@ def mark_boundaries(image, label_img, color=(1, 1, 0),
     if mode == 'subpixel':
         # Here, we want to interpose an extra line of pixels between
         # each original line - except for the last axis which holds
-        # the RGB information. ``nd.zoom`` then performs the (cubic)
+        # the RGB information. ``ndi.zoom`` then performs the (cubic)
         # interpolation, filling in the values of the interposed pixels
-        marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1],
-                         mode='reflect')
+        marked = ndi.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1],
+                          mode='reflect')
     boundaries = find_boundaries(label_img, mode=mode,
                                  background=background_label)
     if outline_color is not None:
diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py
index be325fc6..d3a6a759 100644
--- a/skimage/segmentation/random_walker_segmentation.py
+++ b/skimage/segmentation/random_walker_segmentation.py
@@ -10,7 +10,7 @@ significantly the performance.
 
 import warnings
 import numpy as np
-from scipy import sparse, ndimage
+from scipy import sparse, ndimage as ndi
 
 # executive summary for next code block: try to import umfpack from
 # scipy, but make sure not to raise a fuss if it fails since it's only
@@ -417,7 +417,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True,
     # If the array has pruned zones, be sure that no isolated pixels
     # exist between pruned zones (they could not be determined)
     if np.any(labels < 0):
-        filled = ndimage.binary_propagation(labels > 0, mask=labels >= 0)
+        filled = ndi.binary_propagation(labels > 0, mask=labels >= 0)
         labels[np.logical_and(np.logical_not(filled), labels == 0)] = -1
         del filled
     labels = np.atleast_3d(labels)
diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 08f4f828..ed06272a 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -2,7 +2,7 @@
 
 import collections as coll
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 import warnings
 
 from ..util import img_as_float, regular_grid
@@ -139,7 +139,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0,
     if (sigma > 0).any():
         # add zero smoothing for multichannel dimension
         sigma = list(sigma) + [0]
-        image = ndimage.gaussian_filter(image, sigma)
+        image = ndi.gaussian_filter(image, sigma)
 
     if multichannel and (convert2lab or convert2lab is None):
         if image.shape[-1] != 3 and convert2lab:
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 33448353..073066dd 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -2,11 +2,10 @@ import six
 import math
 import warnings
 import numpy as np
-from scipy import ndimage, spatial
+from scipy import ndimage as ndi, spatial
 
 from .._shared.utils import get_bound_method_class, safe_as_int
 from ..util import img_as_float
-from ..exposure import rescale_intensity
 from ._warps_cy import _warp_fast
 
 
@@ -1054,7 +1053,7 @@ def warp_coords(coord_map, shape, dtype=np.float64):
     users who would like, for example, to re-use a particular coordinate
     mapping, to use specific dtypes at various points along the the
     image-warping process, or to implement different post-processing logic
-    than `warp` performs after the call to `ndimage.map_coordinates`.
+    than `warp` performs after the call to `ndi.map_coordinates`.
 
 
     Examples
@@ -1352,7 +1351,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
                 warped = np.dstack(dims)
 
     if warped is None:
-        # use ndimage.map_coordinates
+        # use ndi.map_coordinates
 
         if (isinstance(inverse_map, np.ndarray)
                 and inverse_map.shape == (3, 3)):
@@ -1388,8 +1387,8 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
         # Pre-filtering not necessary for order 0, 1 interpolation
         prefilter = order > 1
 
-        warped = ndimage.map_coordinates(image, coords, prefilter=prefilter,
-                                      mode=mode, order=order, cval=cval)
+        warped = ndi.map_coordinates(image, coords, prefilter=prefilter,
+                                     mode=mode, order=order, cval=cval)
 
 
     _clip_warp_output(image, warped, order, mode, cval, clip)
diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py
index 5a9c8aa2..f069122a 100644
--- a/skimage/transform/_warps.py
+++ b/skimage/transform/_warps.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 
 from ..measure import block_reduce
 from ._geometric import (warp, SimilarityTransform, AffineTransform,
@@ -81,8 +81,8 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True,
 
         image = _convert_warp_input(image, preserve_range)
 
-        out = ndimage.map_coordinates(image, coord_map, order=order,
-                                      mode=mode, cval=cval)
+        out = ndi.map_coordinates(image, coord_map, order=order,
+                                  mode=mode, cval=cval)
 
         _clip_warp_output(image, out, order, mode, cval, clip)
 
diff --git a/skimage/transform/hough_transform.py b/skimage/transform/hough_transform.py
index 0760f5e7..363caf2e 100644
--- a/skimage/transform/hough_transform.py
+++ b/skimage/transform/hough_transform.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from .. import measure, morphology
 from ._hough_transform import _hough_circle
 
@@ -63,10 +63,10 @@ def hough_line_peaks(hspace, angles, dists, min_distance=9, min_angle=10,
 
     distance_size = 2 * min_distance + 1
     angle_size = 2 * min_angle + 1
-    hspace_max = ndimage.maximum_filter1d(hspace, size=distance_size, axis=0,
-                                          mode='constant', cval=0)
-    hspace_max = ndimage.maximum_filter1d(hspace_max, size=angle_size, axis=1,
-                                          mode='constant', cval=0)
+    hspace_max = ndi.maximum_filter1d(hspace, size=distance_size, axis=0,
+                                      mode='constant', cval=0)
+    hspace_max = ndi.maximum_filter1d(hspace_max, size=angle_size, axis=1,
+                                      mode='constant', cval=0)
     mask = (hspace == hspace_max)
     hspace *= mask
     hspace_t = hspace > threshold
diff --git a/skimage/transform/pyramids.py b/skimage/transform/pyramids.py
index 435d60b9..fe2d26d1 100644
--- a/skimage/transform/pyramids.py
+++ b/skimage/transform/pyramids.py
@@ -1,6 +1,6 @@
 import math
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage as ndi
 from ..transform import resize
 from ..util import img_as_float
 
@@ -13,12 +13,12 @@ def _smooth(image, sigma, mode, cval):
     # apply Gaussian filter to all dimensions independently
     if image.ndim == 3:
         for dim in range(image.shape[2]):
-            ndimage.gaussian_filter(image[..., dim], sigma,
-                                    output=smoothed[..., dim],
-                                    mode=mode, cval=cval)
-    else:
-        ndimage.gaussian_filter(image, sigma, output=smoothed,
+            ndi.gaussian_filter(image[..., dim], sigma,
+                                output=smoothed[..., dim],
                                 mode=mode, cval=cval)
+    else:
+        ndi.gaussian_filter(image, sigma, output=smoothed,
+                            mode=mode, cval=cval)
 
     return smoothed