diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py
index 34e05352..5b2b92db 100644
--- a/skimage/feature/greycomatrix.py
+++ b/skimage/feature/greycomatrix.py
@@ -181,11 +181,11 @@ def greycoprops(P, prop='contrast'):
     # create weights for specified property
     I, J = np.ogrid[0:num_level, 0:num_level]
     if prop == 'contrast':
-        weights = (I - J) ** 2
+        weights = (I - J)**2
     elif prop == 'dissimilarity':
         weights = np.abs(I - J)
     elif prop == 'homogeneity':
-        weights = 1. / (1. + (I - J) ** 2)
+        weights = 1. / (1. + (I - J)**2)
     elif prop in ['ASM', 'energy', 'correlation']:
         pass
     else:
@@ -193,10 +193,10 @@ def greycoprops(P, prop='contrast'):
 
     # compute property for each GLCM
     if prop == 'energy':
-        asm = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
+        asm = np.apply_over_axes(np.sum, (P**2), axes=(0, 1))[0, 0]
         results = np.sqrt(asm)
     elif prop == 'ASM':
-        results = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
+        results = np.apply_over_axes(np.sum, (P**2), axes=(0, 1))[0, 0]
     elif prop == 'correlation':
         results = np.zeros((num_dist, num_angle), dtype=np.float64)
         I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
@@ -204,9 +204,9 @@ def greycoprops(P, prop='contrast'):
         diff_i = I - np.apply_over_axes(np.sum, (I * P), axes=(0, 1))[0, 0]
         diff_j = J - np.apply_over_axes(np.sum, (J * P), axes=(0, 1))[0, 0]
 
-        std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i) ** 2),
+        std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i)**2),
                                            axes=(0, 1))[0, 0])
-        std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j) ** 2),
+        std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j)**2),
                                            axes=(0, 1))[0, 0])
         cov = np.apply_over_axes(np.sum, (P * (diff_i * diff_j)),
                                  axes=(0, 1))[0, 0]
diff --git a/skimage/feature/harris.py b/skimage/feature/harris.py
index e496b33d..14854ac4 100644
--- a/skimage/feature/harris.py
+++ b/skimage/feature/harris.py
@@ -42,7 +42,7 @@ def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1):
     Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant')
 
     # determinant and trace
-    Wdet = Wxx * Wyy - Wxy ** 2
+    Wdet = Wxx * Wyy - Wxy**2
     Wtr = Wxx + Wyy
     # Alternate formula for Harris response.
     # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989)
diff --git a/skimage/feature/hog.py b/skimage/feature/hog.py
index e0e1301d..4a24b0a2 100644
--- a/skimage/feature/hog.py
+++ b/skimage/feature/hog.py
@@ -95,7 +95,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
     cell are used to vote into the orientation histogram.
     """
 
-    magnitude = sqrt(gx ** 2 + gy ** 2)
+    magnitude = sqrt(gx**2 + gy**2)
     orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90
 
     sy, sx = image.shape
@@ -166,7 +166,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
         for y in range(n_blocksy):
             block = orientation_histogram[y:y + by, x:x + bx, :]
             eps = 1e-5
-            normalised_blocks[y, x, :] = block / sqrt(block.sum() ** 2 + eps)
+            normalised_blocks[y, x, :] = block / sqrt(block.sum()**2 + eps)
 
     """
     The final step collects the HOG descriptors from all blocks of a dense
diff --git a/skimage/filter/edges.py b/skimage/filter/edges.py
index 6ceb5c5f..b6a800c9 100644
--- a/skimage/filter/edges.py
+++ b/skimage/filter/edges.py
@@ -37,7 +37,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.
     """
-    return np.sqrt(hsobel(image, mask) ** 2 + vsobel(image, mask) ** 2)
+    return np.sqrt(hsobel(image, mask)**2 + vsobel(image, mask)**2)
 
 
 def hsobel(image, mask=None):
@@ -137,7 +137,7 @@ def prewitt(image, mask=None):
     Return the square root of the sum of squares of the horizontal
     and vertical Prewitt transforms.
     """
-    return np.sqrt(hprewitt(image, mask) ** 2 + vprewitt(image, mask) ** 2)
+    return np.sqrt(hprewitt(image, mask)**2 + vprewitt(image, mask)**2)
 
 
 def hprewitt(image, mask=None):
diff --git a/skimage/filter/lpi_filter.py b/skimage/filter/lpi_filter.py
index 92490abe..60eb1d63 100644
--- a/skimage/filter/lpi_filter.py
+++ b/skimage/filter/lpi_filter.py
@@ -232,7 +232,7 @@ def wiener(data, impulse_response=None, filter_params={}, K=0.25,
     F, G = filt._prepare(data)
     _min_limit(F)
 
-    H_mag_sqr = np.abs(F) ** 2
+    H_mag_sqr = np.abs(F)**2
     F = 1 / F * H_mag_sqr / (H_mag_sqr + K)
 
     return _centre(np.abs(ifftshift(np.dual.ifftn(G * F))), data.shape)
diff --git a/skimage/filter/tests/test_tv_denoise.py b/skimage/filter/tests/test_tv_denoise.py
index 27db4894..4cc6adbb 100644
--- a/skimage/filter/tests/test_tv_denoise.py
+++ b/skimage/filter/tests/test_tv_denoise.py
@@ -27,8 +27,8 @@ class TestTvDenoise():
         grad_denoised = ndimage.morphological_gradient(
             denoised_lena, size=((3, 3)))
         # test if the total variation has decreased
-        assert (np.sqrt((grad_denoised ** 2).sum())
-                < np.sqrt((grad ** 2).sum()) / 2)
+        assert (np.sqrt((grad_denoised**2).sum())
+                < np.sqrt((grad**2).sum()) / 2)
         denoised_lena_int = filter.tv_denoise(img_as_uint(lena),
                                               weight=60.0, keep_type=True)
         assert denoised_lena_int.dtype is np.dtype('uint16')
@@ -39,7 +39,7 @@ class TestTvDenoise():
         a sphere.
         """
         x, y, z = np.ogrid[0:40, 0:40, 0:40]
-        mask = (x - 22) ** 2 + (y - 20) ** 2 + (z - 17) ** 2 < 8 ** 2
+        mask = (x - 22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
         mask = 100 * mask.astype(np.float)
         mask += 60
         mask += 20 * np.random.randn(*mask.shape)
diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py
index 505123b7..022cacc5 100644
--- a/skimage/filter/thresholding.py
+++ b/skimage/filter/thresholding.py
@@ -127,7 +127,7 @@ def threshold_otsu(image, nbins=256):
     # Clip ends to align class 1 and class 2 variables:
     # The last value of `weight1`/`mean1` should pair with zero values in
     # `weight2`/`mean2`, which do not exist.
-    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2
+    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:])**2
 
     idx = np.argmax(variance12)
     threshold = bin_centers[:-1][idx]
diff --git a/skimage/filter/tv_denoise.py b/skimage/filter/tv_denoise.py
index 37abdbb5..c8e3216c 100644
--- a/skimage/filter/tv_denoise.py
+++ b/skimage/filter/tv_denoise.py
@@ -56,12 +56,12 @@ def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
         d[:, :, 1:] += pz[:, :, :-1]
 
         out = im + d
-        E = (d ** 2).sum()
+        E = (d**2).sum()
 
         gx[:-1] = np.diff(out, axis=0)
         gy[:, :-1] = np.diff(out, axis=1)
         gz[:, :, :-1] = np.diff(out, axis=2)
-        norm = np.sqrt(gx ** 2 + gy ** 2 + gz ** 2)
+        norm = np.sqrt(gx**2 + gy**2 + gz**2)
         E += weight * norm.sum()
         norm *= 0.5 / weight
         norm += 1.
@@ -147,10 +147,10 @@ def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
         d[:, 1:] += py[:, :-1]
 
         out = im + d
-        E = (d ** 2).sum()
+        E = (d**2).sum()
         gx[:-1] = np.diff(out, axis=0)
         gy[:, :-1] = np.diff(out, axis=1)
-        norm = np.sqrt(gx ** 2 + gy ** 2)
+        norm = np.sqrt(gx**2 + gy**2)
         E += weight * norm.sum()
         norm *= 0.5 / weight
         norm += 1
diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py
index 89dbddce..6b8f6a3c 100644
--- a/skimage/measure/_regionprops.py
+++ b/skimage/measure/_regionprops.py
@@ -210,8 +210,8 @@ def regionprops(label_image, properties=['Area', 'Centroid'],
         b = mu[1, 1] / mu[0, 0]
         c = mu[0, 2] / mu[0, 0]
         #: eigen values of inertia tensor
-        l1 = (a + c) / 2 + sqrt(4 * b ** 2 + (a - c) ** 2) / 2
-        l2 = (a + c) / 2 - sqrt(4 * b ** 2 + (a - c) ** 2) / 2
+        l1 = (a + c) / 2 + sqrt(4 * b**2 + (a - c)**2) / 2
+        l2 = (a + c) / 2 - sqrt(4 * b**2 + (a - c)**2) / 2
 
         # cached results which are used by several properties
         _filled_image = None
diff --git a/skimage/measure/tests/test_find_contours.py b/skimage/measure/tests/test_find_contours.py
index 1c6096ee..3de7acc4 100644
--- a/skimage/measure/tests/test_find_contours.py
+++ b/skimage/measure/tests/test_find_contours.py
@@ -17,7 +17,7 @@ a[1, 1:-1] = 0
 ##        [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]], dtype=float32)
 
 x, y = np.mgrid[-1:1:5j, -1:1:5j]
-r = np.sqrt(x ** 2 + y ** 2)
+r = np.sqrt(x**2 + y**2)
 
 
 def test_binary():
diff --git a/skimage/morphology/selem.py b/skimage/morphology/selem.py
index 5ff6793a..f2234ffb 100644
--- a/skimage/morphology/selem.py
+++ b/skimage/morphology/selem.py
@@ -111,6 +111,6 @@ def disk(radius, dtype=np.uint8):
     """
     L = np.linspace(-radius, radius, 2 * radius + 1)
     (X, Y) = np.meshgrid(L, L)
-    s = X ** 2
-    s += Y ** 2
+    s = X**2
+    s += Y**2
     return np.array(s <= radius * radius, dtype=dtype)
diff --git a/skimage/morphology/skeletonize.py b/skimage/morphology/skeletonize.py
index b7c59ac1..bfa7525c 100644
--- a/skimage/morphology/skeletonize.py
+++ b/skimage/morphology/skeletonize.py
@@ -244,11 +244,11 @@ def medial_axis(image, mask=None, return_distance=False):
     # OR
     # 3. Keep if # pixels in neighbourhood is 2 or less
     # Note that table is independent of image
-    center_is_foreground = (np.arange(512) & 2 ** 4).astype(bool)
+    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),
+                        ndimage.label(_pattern_of(index & ~ 2**4),
                                     _eight_connect)[1]
                         for index in range(512)])  # condition 2
                 |
@@ -311,9 +311,9 @@ def _pattern_of(index):
     Return the pattern represented by an index value
     Byte decomposition of index
     """
-    return np.array([[index & 2 ** 0, index & 2 ** 1, index & 2 ** 2],
-                     [index & 2 ** 3, index & 2 ** 4, index & 2 ** 5],
-                     [index & 2 ** 6, index & 2 ** 7, index & 2 ** 8]], bool)
+    return np.array([[index & 2**0, index & 2**1, index & 2**2],
+                     [index & 2**3, index & 2**4, index & 2**5],
+                     [index & 2**6, index & 2**7, index & 2**8]], bool)
 
 
 def _table_lookup(image, table):
diff --git a/skimage/morphology/tests/test_skeletonize.py b/skimage/morphology/tests/test_skeletonize.py
index 408e971c..2f9d046e 100644
--- a/skimage/morphology/tests/test_skeletonize.py
+++ b/skimage/morphology/tests/test_skeletonize.py
@@ -80,8 +80,8 @@ class TestSkeletonize():
 
         # foreground object 3
         ir, ic = np.indices(image.shape)
-        circle1 = (ic - 135) ** 2 + (ir - 150) ** 2 < 30 ** 2
-        circle2 = (ic - 135) ** 2 + (ir - 150) ** 2 < 20 ** 2
+        circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2
+        circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
         image[circle1] = 1
         image[circle2] = 0
         result = skeletonize(image)
diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py
index 1fe8baba..c6671d7a 100644
--- a/skimage/morphology/tests/test_watershed.py
+++ b/skimage/morphology/tests/test_watershed.py
@@ -72,7 +72,7 @@ def diff(a, b):
         a = a.astype(np.float64)
         b = np.asarray(b)
         b = b.astype(np.float64)
-        t = ((a - b) ** 2).sum()
+        t = ((a - b)**2).sum()
     return math.sqrt(t)
 
 
diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py
index 130be926..68afcc65 100644
--- a/skimage/segmentation/random_walker_segmentation.py
+++ b/skimage/segmentation/random_walker_segmentation.py
@@ -63,7 +63,7 @@ def _make_graph_edges_3d(n_x, n_y, n_z):
 
 
 def _compute_weights_3d(data, beta=130, eps=1.e-6):
-    gradients = _compute_gradients_3d(data) ** 2
+    gradients = _compute_gradients_3d(data)**2
     beta /= 10 * data.std()
     gradients *= beta
     weights = np.exp(- gradients)
diff --git a/skimage/transform/_warp_zoo.py b/skimage/transform/_warp_zoo.py
index c2e8b4df..4df531d1 100644
--- a/skimage/transform/_warp_zoo.py
+++ b/skimage/transform/_warp_zoo.py
@@ -7,7 +7,7 @@ from ._warp import warp
 def _swirl_mapping(xy, center, rotation, strength, radius):
     x, y = xy.T
     x0, y0 = center
-    rho = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
+    rho = np.sqrt((x - x0)**2 + (y - y0)**2)
 
     # Ensure that the transformation decays to approximately 1/1000-th
     # within the specified radius.
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 2480c7fc..d185c0fc 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -43,7 +43,7 @@ def radon(image, theta=None):
     if theta == None:
         theta = np.arange(180)
     height, width = image.shape
-    diagonal = np.sqrt(height ** 2 + width ** 2)
+    diagonal = np.sqrt(height**2 + width**2)
     heightpad = np.ceil(diagonal - height)
     widthpad = np.ceil(diagonal - width)
     padded_image = np.zeros((int(height + heightpad),
@@ -130,13 +130,13 @@ def iradon(radon_image, theta=None, output_size=None,
     th = (np.pi / 180.0) * theta
     # if output size not specified, estimate from input radon image
     if not output_size:
-        output_size = int(np.floor(np.sqrt((radon_image.shape[0]) ** 2 / 2.0)))
+        output_size = int(np.floor(np.sqrt((radon_image.shape[0])**2 / 2.0)))
     n = radon_image.shape[0]
 
     img = radon_image.copy()
     # resize image to next power of two for fourier analysis
     # speeds up fourier and lessens artifacts
-    order = max(64., 2 ** np.ceil(np.log(2 * n) / np.log(2)))
+    order = max(64., 2**np.ceil(np.log(2 * n) / np.log(2)))
     # zero pad input image
     img.resize((order, img.shape[1]))
     # construct the fourier filter
diff --git a/skimage/util/dtype.py b/skimage/util/dtype.py
index 6383c53f..9697e627 100644
--- a/skimage/util/dtype.py
+++ b/skimage/util/dtype.py
@@ -109,21 +109,21 @@ def convert(image, dtype, force_copy=False, uniform=False):
             prec_loss()
             if copy:
                 b = np.empty(a.shape, _dtype2(kind, m))
-                np.floor_divide(a, 2 ** (n - m), out=b, dtype=a.dtype,
+                np.floor_divide(a, 2**(n - m), out=b, dtype=a.dtype,
                                 casting='unsafe')
                 return b
             else:
-                a //= 2 ** (n - m)
+                a //= 2**(n - m)
                 return a
         elif m % n == 0:
             # exact upscale to a multiple of n bits
             if copy:
                 b = np.empty(a.shape, _dtype2(kind, m))
-                np.multiply(a, (2 ** m - 1) // (2 ** n - 1), out=b, dtype=b.dtype)
+                np.multiply(a, (2**m - 1) // (2**n - 1), out=b, dtype=b.dtype)
                 return b
             else:
                 a = np.array(a, _dtype2(kind, m, a.dtype.itemsize), copy=False)
-                a *= (2 ** m - 1) // (2 ** n - 1)
+                a *= (2**m - 1) // (2**n - 1)
                 return a
         else:
             # upscale to a multiple of n bits,
@@ -132,13 +132,13 @@ def convert(image, dtype, force_copy=False, uniform=False):
             o = (m // n + 1) * n
             if copy:
                 b = np.empty(a.shape, _dtype2(kind, o))
-                np.multiply(a, (2 ** o - 1) // (2 ** n - 1), out=b, dtype=b.dtype)
-                b //= 2 ** (o - m)
+                np.multiply(a, (2**o - 1) // (2**n - 1), out=b, dtype=b.dtype)
+                b //= 2**(o - m)
                 return b
             else:
                 a = np.array(a, _dtype2(kind, o, a.dtype.itemsize), copy=False)
-                a *= (2 ** o - 1) // (2 ** n - 1)
-                a //= 2 ** (o - m)
+                a *= (2**o - 1) // (2**n - 1)
+                a //= 2**(o - m)
                 return a
 
     kind = dtypeobj.kind
diff --git a/skimage/util/montage.py b/skimage/util/montage.py
index 57ec7274..6c23ccbf 100644
--- a/skimage/util/montage.py
+++ b/skimage/util/montage.py
@@ -84,7 +84,7 @@ def montage2d(arr_in, fill='mean', rescale_intensity=False):
     if fill == 'mean':
         fill = arr_in.mean()
 
-    n_missing = int((alpha ** 2.) - n_images)
+    n_missing = int((alpha**2.) - n_images)
     missing = np.ones((n_missing, height, width), dtype=arr_in.dtype) * fill
     arr_out = np.vstack((arr_in, missing))