From 5c09f39789ce50aa4bdff55157768484d299a900 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 29 Apr 2021 14:52:14 -0400
Subject: [PATCH 1/8] float32 support in mark_boundaries

---
 skimage/segmentation/boundaries.py            |  3 +++
 skimage/segmentation/tests/test_boundaries.py | 13 ++++++++++---
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py
index 562387c2ee5..da466eec4ff 100644
--- a/skimage/segmentation/boundaries.py
+++ b/skimage/segmentation/boundaries.py
@@ -1,6 +1,7 @@
 
 import numpy as np
 from scipy import ndimage as ndi
+from .._shared.utils import _supported_float_type
 from ..morphology import dilation, erosion, square
 from ..util import img_as_float, view_as_windows
 from ..color import gray2rgb
@@ -212,7 +213,9 @@ def mark_boundaries(image, label_img, color=(1, 1, 0),
     --------
     find_boundaries
     """
+    float_dtype = _supported_float_type(image)
     marked = img_as_float(image, force_copy=True)
+    marked = marked.astype(float_dtype, copy=False)
     if marked.ndim == 2:
         marked = gray2rgb(marked)
     if mode == 'subpixel':
diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py
index 3c6577219c2..e291831c15f 100644
--- a/skimage/segmentation/tests/test_boundaries.py
+++ b/skimage/segmentation/tests/test_boundaries.py
@@ -1,7 +1,9 @@
 import numpy as np
 from skimage.segmentation import find_boundaries, mark_boundaries
 
+from skimage._shared import testing
 from skimage._shared.testing import assert_array_equal, assert_allclose
+from skimage._shared.utils import _supported_float_type
 
 
 white = (1, 1, 1)
@@ -39,8 +41,9 @@ def test_find_boundaries_bool():
     assert_array_equal(result, ref)
 
 
-def test_mark_boundaries():
-    image = np.zeros((10, 10))
+@testing.parametrize('dtype', [np.uint8, np.float16, np.float32, np.float64])
+def test_mark_boundaries(dtype):
+    image = np.zeros((10, 10), dtype=dtype)
     label_image = np.zeros((10, 10), dtype=np.uint8)
     label_image[2:7, 2:7] = 1
 
@@ -56,6 +59,7 @@ def test_mark_boundaries():
                     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
 
     marked = mark_boundaries(image, label_image, color=white, mode='thick')
+    assert marked.dtype == _supported_float_type(dtype)
     result = np.mean(marked, axis=-1)
     assert_array_equal(result, ref)
 
@@ -96,7 +100,8 @@ def test_mark_boundaries_bool():
     assert_array_equal(result, ref)
 
 
-def test_mark_boundaries_subpixel():
+@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_mark_boundaries_subpixel(dtype):
     labels = np.array([[0, 0, 0, 0],
                        [0, 0, 5, 0],
                        [0, 1, 5, 0],
@@ -104,7 +109,9 @@ def test_mark_boundaries_subpixel():
                        [0, 0, 0, 0]], dtype=np.uint8)
     np.random.seed(0)
     image = np.round(np.random.rand(*labels.shape), 2)
+    image = image.astype(dtype, copy=False)
     marked = mark_boundaries(image, labels, color=white, mode='subpixel')
+    assert marked.dtype == _supported_float_type(dtype)
     marked_proj = np.round(np.mean(marked, axis=-1), 2)
 
     ref_result = np.array(

From e10b77b985f10f6d796d81611de0f95516e1e014 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 29 Apr 2021 15:27:59 -0400
Subject: [PATCH 2/8] float32 support in active_contour

---
 skimage/segmentation/active_contour_model.py  | 54 +++++++++++--------
 skimage/segmentation/boundaries.py            |  1 +
 .../tests/test_active_contour_model.py        | 35 ++++++++----
 3 files changed, 59 insertions(+), 31 deletions(-)

diff --git a/skimage/segmentation/active_contour_model.py b/skimage/segmentation/active_contour_model.py
index 18621cbfb0a..7f4d1e71224 100644
--- a/skimage/segmentation/active_contour_model.py
+++ b/skimage/segmentation/active_contour_model.py
@@ -1,5 +1,7 @@
 import numpy as np
 from scipy.interpolate import RectBivariateSpline
+
+from .._shared.utils import _supported_float_type
 from ..util import img_as_float
 from ..filters import sobel
 
@@ -108,7 +110,11 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
     if boundary_condition not in valid_bcs:
         raise ValueError("Invalid boundary condition.\n" +
                          "Should be one of: "+", ".join(valid_bcs)+'.')
+
     img = img_as_float(image)
+    float_dtype = _supported_float_type(image)
+    img = img.astype(float_dtype, copy=False)
+
     RGB = img.ndim == 3
 
     # Find edges using sobel:
@@ -134,20 +140,22 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
                                img.T, kx=2, ky=2, s=0)
 
     snake_xy = snake[:, ::-1]
-    x, y = snake_xy[:, 0].astype(float), snake_xy[:, 1].astype(float)
+    x = snake_xy[:, 0].astype(float_dtype)
+    y = snake_xy[:, 1].astype(float_dtype)
     n = len(x)
-    xsave = np.empty((convergence_order, n))
-    ysave = np.empty((convergence_order, n))
-
-    # Build snake shape matrix for Euler equation
-    a = np.roll(np.eye(n), -1, axis=0) + \
-        np.roll(np.eye(n), -1, axis=1) - \
-        2*np.eye(n)  # second order derivative, central difference
-    b = np.roll(np.eye(n), -2, axis=0) + \
-        np.roll(np.eye(n), -2, axis=1) - \
-        4*np.roll(np.eye(n), -1, axis=0) - \
-        4*np.roll(np.eye(n), -1, axis=1) + \
-        6*np.eye(n)  # fourth order derivative, central difference
+    xsave = np.empty((convergence_order, n), dtype=float_dtype)
+    ysave = np.empty((convergence_order, n), dtype=float_dtype)
+
+    # Build snake shape matrix for Euler equation in double precision
+    eye_n = np.eye(n, dtype=float)
+    a = np.roll(eye_n, -1, axis=0) + \
+        np.roll(eye_n, -1, axis=1) - \
+        2*eye_n  # second order derivative, central difference
+    b = np.roll(eye_n, -2, axis=0) + \
+        np.roll(eye_n, -2, axis=1) - \
+        4*np.roll(eye_n, -1, axis=0) - \
+        4*np.roll(eye_n, -1, axis=1) + \
+        6*eye_n  # fourth order derivative, central difference
     A = -alpha*a + beta*b
 
     # Impose boundary conditions different from periodic:
@@ -179,12 +187,16 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
         efree = True
 
     # Only one inversion is needed for implicit spline energy minimization:
-    inv = np.linalg.inv(A + gamma*np.eye(n))
+    inv = np.linalg.inv(A + gamma*eye_n)
+    # can use float_dtype once we have computed the inverse in double precision
+    inv = inv.astype(float_dtype, copy=False)
 
     # Explicit time stepping for image energy minimization:
     for i in range(max_iterations):
-        fx = intp(x, y, dx=1, grid=False)
-        fy = intp(x, y, dy=1, grid=False)
+        # RectBivariateSpline always returns float64, so call astype here
+        fx = intp(x, y, dx=1, grid=False).astype(float_dtype, copy=False)
+        fy = intp(x, y, dy=1, grid=False).astype(float_dtype, copy=False)
+
         if sfixed:
             fx[0] = 0
             fy[0] = 0
@@ -201,8 +213,8 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
         yn = inv @ (gamma*y + fy)
 
         # Movements are capped to max_px_move per iteration:
-        dx = max_px_move*np.tanh(xn-x)
-        dy = max_px_move*np.tanh(yn-y)
+        dx = max_px_move*np.tanh(xn - x)
+        dy = max_px_move*np.tanh(yn - y)
         if sfixed:
             dx[0] = 0
             dy[0] = 0
@@ -214,13 +226,13 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
 
         # Convergence criteria needs to compare to a number of previous
         # configurations since oscillations can occur.
-        j = i % (convergence_order+1)
+        j = i % (convergence_order + 1)
         if j < convergence_order:
             xsave[j, :] = x
             ysave[j, :] = y
         else:
-            dist = np.min(np.max(np.abs(xsave-x[None, :]) +
-                                 np.abs(ysave-y[None, :]), 1))
+            dist = np.min(np.max(np.abs(xsave - x[None, :]) +
+                                 np.abs(ysave - y[None, :]), 1))
             if dist < convergence:
                 break
 
diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py
index da466eec4ff..d3937cbe05c 100644
--- a/skimage/segmentation/boundaries.py
+++ b/skimage/segmentation/boundaries.py
@@ -1,6 +1,7 @@
 
 import numpy as np
 from scipy import ndimage as ndi
+
 from .._shared.utils import _supported_float_type
 from ..morphology import dilation, erosion, square
 from ..util import img_as_float, view_as_windows
diff --git a/skimage/segmentation/tests/test_active_contour_model.py b/skimage/segmentation/tests/test_active_contour_model.py
index 1fed16af09e..f4b353cc1b2 100644
--- a/skimage/segmentation/tests/test_active_contour_model.py
+++ b/skimage/segmentation/tests/test_active_contour_model.py
@@ -6,54 +6,65 @@
 
 from skimage._shared import testing
 from skimage._shared.testing import assert_equal, assert_allclose
+from skimage._shared.utils import _supported_float_type
 
 
-def test_periodic_reference():
+@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_periodic_reference(dtype):
     img = data.astronaut()
     img = rgb2gray(img)
     s = np.linspace(0, 2*np.pi, 400)
     r = 100 + 100*np.sin(s)
     c = 220 + 100*np.cos(s)
     init = np.array([r, c]).T
-    snake = active_contour(gaussian(img, 3), init, alpha=0.015, beta=10,
+    img_smooth = gaussian(img, 3).astype(dtype, copy=False)
+    snake = active_contour(img_smooth, init, alpha=0.015, beta=10,
                            w_line=0, w_edge=1, gamma=0.001)
+    assert snake.dtype == _supported_float_type(dtype)
     refr = [98, 99, 100, 101, 102, 103, 104, 105, 106, 108]
     refc = [299, 298, 298, 298, 298, 297, 297, 296, 296, 295]
     assert_equal(np.array(snake[:10, 0], dtype=np.int32), refr)
     assert_equal(np.array(snake[:10, 1], dtype=np.int32), refc)
 
 
-def test_fixed_reference():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_fixed_reference(dtype):
     img = data.text()
     r = np.linspace(136, 50, 100)
     c = np.linspace(5, 424, 100)
     init = np.array([r, c]).T
-    snake = active_contour(gaussian(img, 1), init, boundary_condition='fixed',
+    image_smooth = gaussian(img, 1).astype(dtype, copy=False)
+    snake = active_contour(image_smooth, init, boundary_condition='fixed',
                            alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
+    assert snake.dtype == _supported_float_type(dtype)
     refr = [136, 135, 134, 133, 132, 131, 129, 128, 127, 125]
     refc = [5, 9, 13, 17, 21, 25, 30, 34, 38, 42]
     assert_equal(np.array(snake[:10, 0], dtype=np.int32), refr)
     assert_equal(np.array(snake[:10, 1], dtype=np.int32), refc)
 
 
-def test_free_reference():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_free_reference(dtype):
     img = data.text()
     r = np.linspace(70, 40, 100)
     c = np.linspace(5, 424, 100)
     init = np.array([r, c]).T
-    snake = active_contour(gaussian(img, 3), init, boundary_condition='free',
+    img_smooth = gaussian(img, 3).astype(dtype, copy=False)
+    snake = active_contour(img_smooth, init, boundary_condition='free',
                            alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
+    assert snake.dtype == _supported_float_type(dtype)
     refr = [76, 76, 75, 74, 73, 72, 71, 70, 69, 69]
     refc = [10, 13, 16, 19, 23, 26, 29, 32, 36, 39]
     assert_equal(np.array(snake[:10, 0], dtype=np.int32), refr)
     assert_equal(np.array(snake[:10, 1], dtype=np.int32), refc)
 
 
-def test_RGB():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_RGB(dtype):
     img = gaussian(data.text(), 1)
-    imgR = np.zeros((img.shape[0], img.shape[1], 3))
-    imgG = np.zeros((img.shape[0], img.shape[1], 3))
-    imgRGB = np.zeros((img.shape[0], img.shape[1], 3))
+    imgR = np.zeros((img.shape[0], img.shape[1], 3), dtype=dtype)
+    imgG = np.zeros((img.shape[0], img.shape[1], 3), dtype=dtype)
+    imgRGB = np.zeros((img.shape[0], img.shape[1], 3), dtype=dtype)
     imgR[:, :, 0] = img
     imgG[:, :, 1] = img
     imgRGB[:, :, :] = img[:, :, None]
@@ -62,17 +73,21 @@ def test_RGB():
     init = np.array([r, c]).T
     snake = active_contour(imgR, init, boundary_condition='fixed',
                            alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
+    float_dtype = _supported_float_type(dtype)
+    assert snake.dtype == float_dtype
     refr = [136, 135, 134, 133, 132, 131, 129, 128, 127, 125]
     refc = [5, 9, 13, 17, 21, 25, 30, 34, 38, 42]
     assert_equal(np.array(snake[:10, 0], dtype=np.int32), refr)
     assert_equal(np.array(snake[:10, 1], dtype=np.int32), refc)
     snake = active_contour(imgG, init, boundary_condition='fixed',
                            alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1)
+    assert snake.dtype == float_dtype
     assert_equal(np.array(snake[:10, 0], dtype=np.int32), refr)
     assert_equal(np.array(snake[:10, 1], dtype=np.int32), refc)
     snake = active_contour(imgRGB, init, boundary_condition='fixed',
                            alpha=0.1, beta=1.0, w_line=-5/3., w_edge=0,
                            gamma=0.1)
+    assert snake.dtype == float_dtype
     assert_equal(np.array(snake[:10, 0], dtype=np.int32), refr)
     assert_equal(np.array(snake[:10, 1], dtype=np.int32), refc)
 

From 043cbb4eeefc157b6c1a58e42abd34df60d790d8 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 29 Apr 2021 17:39:53 -0400
Subject: [PATCH 3/8] float32 support for chan_vese

---
 skimage/segmentation/_chan_vese.py           | 28 ++++++++++++--------
 skimage/segmentation/tests/test_chan_vese.py | 24 ++++++++++-------
 2 files changed, 32 insertions(+), 20 deletions(-)

diff --git a/skimage/segmentation/_chan_vese.py b/skimage/segmentation/_chan_vese.py
index a8906f632df..44e8f3d2db6 100644
--- a/skimage/segmentation/_chan_vese.py
+++ b/skimage/segmentation/_chan_vese.py
@@ -1,6 +1,8 @@
 import numpy as np
 from scipy.ndimage import distance_transform_edt as distance
 
+from .._shared.utils import _supported_float_type
+
 
 def _cv_curvature(phi):
     """Returns the 'curvature' of a level set 'phi'.
@@ -113,15 +115,17 @@ def _cv_reset_level_set(phi):
     return phi
 
 
-def _cv_checkerboard(image_size, square_size):
+def _cv_checkerboard(image_size, square_size, dtype=np.float64):
     """Generates a checkerboard level set function.
 
     According to Pascal Getreuer, such a level set function has fast convergence.
     """
-    yv = np.arange(image_size[0]).reshape(image_size[0], 1)
-    xv = np.arange(image_size[1])
-    return (np.sin(np.pi/square_size*yv) *
-            np.sin(np.pi/square_size*xv))
+    yv = np.arange(image_size[0], dtype=dtype).reshape(image_size[0], 1)
+    xv = np.arange(image_size[1], dtype=dtype)
+    sf = np.pi / square_size
+    xv *= sf
+    yv *= sf
+    return np.sin(yv) * np.sin(xv)
 
 
 def _cv_large_disk(image_size):
@@ -134,7 +138,7 @@ def _cv_large_disk(image_size):
     centerX = int((image_size[1]-1) / 2)
     res[centerY, centerX] = 0.
     radius = float(min(centerX, centerY))
-    return (radius-distance(res)) / radius
+    return (radius - distance(res)) / radius
 
 
 def _cv_small_disk(image_size):
@@ -147,15 +151,15 @@ def _cv_small_disk(image_size):
     centerX = int((image_size[1]-1) / 2)
     res[centerY, centerX] = 0.
     radius = float(min(centerX, centerY)) / 2.0
-    return (radius-distance(res)) / (radius*3)
+    return (radius - distance(res)) / (radius*3)
 
 
-def _cv_init_level_set(init_level_set, image_shape):
+def _cv_init_level_set(init_level_set, image_shape, dtype=np.float64):
     """Generates an initial level set function conditional on input arguments.
     """
     if type(init_level_set) == str:
         if init_level_set == 'checkerboard':
-            res = _cv_checkerboard(image_shape, 5)
+            res = _cv_checkerboard(image_shape, 5, dtype)
         elif init_level_set == 'disk':
             res = _cv_large_disk(image_shape)
         elif init_level_set == 'small disk':
@@ -164,7 +168,7 @@ def _cv_init_level_set(init_level_set, image_shape):
             raise ValueError("Incorrect name for starting level set preset.")
     else:
         res = init_level_set
-    return res
+    return res.astype(dtype, copy=False)
 
 
 def chan_vese(image, mu=0.25, lambda1=1.0, lambda2=1.0, tol=1e-3, max_iter=500,
@@ -297,12 +301,14 @@ def chan_vese(image, mu=0.25, lambda1=1.0, lambda2=1.0, tol=1e-3, max_iter=500,
     if len(image.shape) != 2:
         raise ValueError("Input image should be a 2D array.")
 
-    phi = _cv_init_level_set(init_level_set, image.shape)
+    float_dtype = _supported_float_type(image)
+    phi = _cv_init_level_set(init_level_set, image.shape, dtype=float_dtype)
 
     if type(phi) != np.ndarray or phi.shape != image.shape:
         raise ValueError("The dimensions of initial level set do not "
                          "match the dimensions of image.")
 
+    image = image.astype(float_dtype, copy=False)
     image = image - np.min(image)
     if np.max(image) != 0:
         image = image / np.max(image)
diff --git a/skimage/segmentation/tests/test_chan_vese.py b/skimage/segmentation/tests/test_chan_vese.py
index 34486967568..2330bcbb7d6 100644
--- a/skimage/segmentation/tests/test_chan_vese.py
+++ b/skimage/segmentation/tests/test_chan_vese.py
@@ -3,18 +3,20 @@
 
 from skimage._shared import testing
 from skimage._shared.testing import assert_array_equal
+from skimage._shared.utils import _supported_float_type
 
 
-def test_chan_vese_flat_level_set():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_chan_vese_flat_level_set(dtype):
     # because the algorithm evolves the level set around the
     # zero-level, it the level-set has no zero level, the algorithm
     # will not produce results in theory. However, since a continuous
     # approximation of the delta function is used, the algorithm
     # still affects the entirety of the level-set. Therefore with
     # infinite time, the segmentation will still converge.
-    img = np.zeros((10, 10))
-    img[3:6, 3:6] = np.ones((3, 3))
-    ls = np.full((10, 10), 1000)
+    img = np.zeros((10, 10), dtype=dtype)
+    img[3:6, 3:6] = 1
+    ls = np.full((10, 10), 1000, dtype=dtype)
     result = chan_vese(img, mu=0.0, tol=1e-3, init_level_set=ls)
     assert_array_equal(result.astype(float), np.ones((10, 10)))
     result = chan_vese(img, mu=0.0, tol=1e-3, init_level_set=-ls)
@@ -23,22 +25,26 @@ def test_chan_vese_flat_level_set():
 
 def test_chan_vese_small_disk_level_set():
     img = np.zeros((10, 10))
-    img[3:6, 3:6] = np.ones((3, 3))
+    img[3:6, 3:6] = 1
     result = chan_vese(img, mu=0.0, tol=1e-3, init_level_set="small disk")
     assert_array_equal(result.astype(float), img)
 
 
 def test_chan_vese_simple_shape():
     img = np.zeros((10, 10))
-    img[3:6, 3:6] = np.ones((3, 3))
+    img[3:6, 3:6] = 1
     result = chan_vese(img, mu=0.0, tol=1e-8).astype(float)
     assert_array_equal(result, img)
 
 
-def test_chan_vese_extended_output():
-    img = np.zeros((10, 10))
-    img[3:6, 3:6] = np.ones((3, 3))
+@testing.parametrize('dtype', [np.uint8, np.float16, np.float32, np.float64])
+def test_chan_vese_extended_output(dtype):
+    img = np.zeros((10, 10), dtype=dtype)
+    img[3:6, 3:6] = 1
     result = chan_vese(img, mu=0.0, tol=1e-8, extended_output=True)
+    float_dtype = _supported_float_type(dtype)
+    assert result[1].dtype == float_dtype
+    assert all(arr.dtype == float_dtype for arr in result[2])
     assert_array_equal(len(result), 3)
 
 

From 538961e36a6d139b6c4e6b853d83f59bd3d5f216 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 29 Apr 2021 18:03:37 -0400
Subject: [PATCH 4/8] float32 support in quickshift

---
 skimage/segmentation/_quickshift.py           |  5 +++-
 skimage/segmentation/_quickshift_cy.pyx       | 26 +++++++++++++------
 skimage/segmentation/tests/test_quickshift.py | 10 +++++--
 3 files changed, 30 insertions(+), 11 deletions(-)

diff --git a/skimage/segmentation/_quickshift.py b/skimage/segmentation/_quickshift.py
index 726618e0e3a..a85063e5940 100644
--- a/skimage/segmentation/_quickshift.py
+++ b/skimage/segmentation/_quickshift.py
@@ -2,9 +2,9 @@
 
 import scipy.ndimage as ndi
 
+from .._shared.utils import _supported_float_type
 from ..util import img_as_float
 from ..color import rgb2lab
-
 from ._quickshift_cy import _quickshift_cython
 
 
@@ -57,6 +57,9 @@ def quickshift(image, ratio=1.0, kernel_size=5, max_dist=10,
     """
 
     image = img_as_float(np.atleast_3d(image))
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
     if convert2lab:
         if image.shape[2] != 3:
             ValueError("Only RGB images can be converted to Lab space.")
diff --git a/skimage/segmentation/_quickshift_cy.pyx b/skimage/segmentation/_quickshift_cy.pyx
index 2abe3262c1e..1be11e32136 100644
--- a/skimage/segmentation/_quickshift_cy.pyx
+++ b/skimage/segmentation/_quickshift_cy.pyx
@@ -6,14 +6,16 @@
 import numpy as np
 cimport numpy as cnp
 
+from .._shared.fused_numerics cimport np_floats
+
 from libc.math cimport exp, sqrt, ceil
 from libc.float cimport DBL_MAX
 
 cnp.import_array()
 
 
-def _quickshift_cython(double[:, :, ::1] image, double kernel_size,
-                       double max_dist, bint return_tree, int random_seed):
+def _quickshift_cython(np_floats[:, :, ::1] image, np_floats kernel_size,
+                       np_floats max_dist, bint return_tree, int random_seed):
     """Segments image using quickshift clustering in Color-(x,y) space.
 
     Produces an oversegmentation of the image using the quickshift mode-seeking
@@ -42,6 +44,11 @@ def _quickshift_cython(double[:, :, ::1] image, double kernel_size,
 
     random_state = np.random.RandomState(random_seed)
 
+    if np_floats is cnp.float64_t:
+        dtype = np.float64
+    else:
+        dtype = np.float32
+
     # TODO join orphaned roots?
     # Some nodes might not have a point of higher density within the
     # search window. We could do a global search over these in the end.
@@ -49,25 +56,28 @@ def _quickshift_cython(double[:, :, ::1] image, double kernel_size,
     # an effect for very high max_dist.
 
     # window size for neighboring pixels to consider
-    cdef double inv_kernel_size_sqr = -0.5 / (kernel_size * kernel_size)
+    cdef np_floats inv_kernel_size_sqr = -0.5 / (kernel_size * kernel_size)
     cdef int kernel_width = <int>ceil(3 * kernel_size)
 
     cdef Py_ssize_t height = image.shape[0]
     cdef Py_ssize_t width = image.shape[1]
     cdef Py_ssize_t channels = image.shape[2]
 
-    cdef double[:, ::1] densities = np.zeros((height, width), dtype=np.double)
+    cdef np_floats[:, ::1] densities = np.zeros((height, width), dtype=dtype)
 
-    cdef double current_density, closest, dist, t
+    cdef np_floats current_density, closest, dist, t
     cdef Py_ssize_t r, c, r_, c_, channel, r_min, r_max, c_min, c_max
-    cdef double* current_pixel_ptr
+    cdef np_floats* current_pixel_ptr
 
     # this will break ties that otherwise would give us headache
-    densities += random_state.normal(scale=0.00001, size=(height, width))
+    densities += random_state.normal(
+        scale=0.00001, size=(height, width)
+    ).astype(dtype, copy=False)
+
     # default parent to self
     cdef Py_ssize_t[:, ::1] parent = \
         np.arange(width * height, dtype=np.intp).reshape(height, width)
-    cdef double[:, ::1] dist_parent = np.zeros((height, width), dtype=np.double)
+    cdef np_floats[:, ::1] dist_parent = np.zeros((height, width), dtype=dtype)
 
     # compute densities
     with nogil:
diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py
index ee0f051b5a0..6c881f3b955 100644
--- a/skimage/segmentation/tests/test_quickshift.py
+++ b/skimage/segmentation/tests/test_quickshift.py
@@ -1,18 +1,22 @@
 import numpy as np
 from skimage.segmentation import quickshift
 
+from skimage._shared import testing
 from skimage._shared.testing import (assert_greater, test_parallel,
                                      assert_equal, assert_array_equal)
+from skimage._shared.utils import _supported_float_type
 
 
 @test_parallel()
-def test_grey():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_grey(dtype):
     rnd = np.random.RandomState(0)
     img = np.zeros((20, 21))
     img[:10, 10:] = 0.2
     img[10:, :10] = 0.4
     img[10:, 10:] = 0.6
     img += 0.1 * rnd.normal(size=img.shape)
+    img = img.astype(dtype, copy=False)
     seg = quickshift(img, kernel_size=2, max_dist=3, random_seed=0,
                      convert2lab=False, sigma=0)
     # we expect 4 segments:
@@ -23,7 +27,8 @@ def test_grey():
         assert_greater(hist[i], 20)
 
 
-def test_color():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_color(dtype):
     rnd = np.random.RandomState(0)
     img = np.zeros((20, 21, 3))
     img[:10, :10, 0] = 1
@@ -32,6 +37,7 @@ def test_color():
     img += 0.01 * rnd.normal(size=img.shape)
     img[img > 1] = 1
     img[img < 0] = 0
+    img = img.astype(dtype, copy=False)
     seg = quickshift(img, random_seed=0, max_dist=30, kernel_size=10, sigma=0)
     # we expect 4 segments:
     assert_equal(len(np.unique(seg)), 4)

From bc492e2287bcead850e7057a3c21f15e43afa5c4 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 29 Apr 2021 18:24:53 -0400
Subject: [PATCH 5/8] add float32 tests for random_walker

---
 .../segmentation/tests/test_random_walker.py  | 38 ++++++++++++++-----
 1 file changed, 28 insertions(+), 10 deletions(-)

diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py
index 4cbb357baa9..8a392451470 100644
--- a/skimage/segmentation/tests/test_random_walker.py
+++ b/skimage/segmentation/tests/test_random_walker.py
@@ -65,16 +65,22 @@ def make_3d_syntheticdata(lx, ly=None, lz=None):
     return data, seeds
 
 
-def test_2d_bf():
+@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_2d_bf(dtype):
     lx = 70
     ly = 100
+
+    # have to use a smaller beta to avoid warning with lower precision input
+    beta = 90 if dtype == np.float64 else 25
+
     data, labels = make_2d_syntheticdata(lx, ly)
+    data = data.astype(dtype, copy=False)
     with expected_warnings([NUMPY_MATRIX_WARNING]):
-        labels_bf = random_walker(data, labels, beta=90, mode='bf')
+        labels_bf = random_walker(data, labels, beta=beta, mode='bf')
     assert (labels_bf[25:45, 40:60] == 2).all()
     assert data.shape == labels.shape
     with expected_warnings([NUMPY_MATRIX_WARNING]):
-        full_prob_bf = random_walker(data, labels, beta=90, mode='bf',
+        full_prob_bf = random_walker(data, labels, beta=beta, mode='bf',
                                      return_full_prob=True)
     assert (full_prob_bf[1, 25:45, 40:60] >=
             full_prob_bf[0, 25:45, 40:60]).all()
@@ -82,7 +88,7 @@ def test_2d_bf():
     # Now test with more than two labels
     labels[55, 80] = 3
     with expected_warnings([NUMPY_MATRIX_WARNING]):
-        full_prob_bf = random_walker(data, labels, beta=90, mode='bf',
+        full_prob_bf = random_walker(data, labels, beta=beta, mode='bf',
                                      return_full_prob=True)
     assert (full_prob_bf[1, 25:45, 40:60] >=
             full_prob_bf[0, 25:45, 40:60]).all()
@@ -90,10 +96,12 @@ def test_2d_bf():
     assert data.shape == labels.shape
 
 
-def test_2d_cg():
+@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_2d_cg(dtype):
     lx = 70
     ly = 100
     data, labels = make_2d_syntheticdata(lx, ly)
+    data = data.astype(dtype, copy=False)
     with expected_warnings(['"cg" mode' + '|' + SCIPY_RANK_WARNING,
                             NUMPY_MATRIX_WARNING]):
         labels_cg = random_walker(data, labels, beta=90, mode='cg')
@@ -109,10 +117,12 @@ def test_2d_cg():
     return data, labels_cg
 
 
-def test_2d_cg_mg():
+@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_2d_cg_mg(dtype):
     lx = 70
     ly = 100
     data, labels = make_2d_syntheticdata(lx, ly)
+    data = data.astype(dtype, copy=False)
     anticipated_warnings = [
         'scipy.sparse.sparsetools|%s' % PYAMG_OR_SCIPY_WARNING,
         NUMPY_MATRIX_WARNING]
@@ -129,10 +139,12 @@ def test_2d_cg_mg():
     return data, labels_cg_mg
 
 
-def test_2d_cg_j():
+@testing.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_2d_cg_j(dtype):
     lx = 70
     ly = 100
     data, labels = make_2d_syntheticdata(lx, ly)
+    data = data.astype(dtype, copy=False)
     with expected_warnings([NUMPY_MATRIX_WARNING]):
         labels_cg = random_walker(data, labels, beta=90, mode='cg_j')
     assert (labels_cg[25:45, 40:60] == 2).all()
@@ -201,10 +213,12 @@ def test_2d_laplacian_size():
     return data, labels
 
 
-def test_3d():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_3d(dtype):
     n = 30
     lx, ly, lz = n, n, n
     data, labels = make_3d_syntheticdata(lx, ly, lz)
+    data = data.astype(dtype, copy=False)
     with expected_warnings(['"cg" mode' + '|' + SCIPY_RANK_WARNING,
                             NUMPY_MATRIX_WARNING]):
         labels = random_walker(data, labels, mode='cg')
@@ -228,9 +242,11 @@ def test_3d_inactive():
     return data, labels, old_labels, after_labels
 
 
-def test_multispectral_2d():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_multispectral_2d(dtype):
     lx, ly = 70, 100
     data, labels = make_2d_syntheticdata(lx, ly)
+    data = data.astype(dtype, copy=False)
     data = data[..., np.newaxis].repeat(2, axis=-1)  # Expect identical output
     with expected_warnings(['"cg" mode' + '|' + SCIPY_RANK_WARNING,
                             NUMPY_MATRIX_WARNING,
@@ -246,10 +262,12 @@ def test_multispectral_2d():
     return data, multi_labels, single_labels, labels
 
 
-def test_multispectral_3d():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_multispectral_3d(dtype):
     n = 30
     lx, ly, lz = n, n, n
     data, labels = make_3d_syntheticdata(lx, ly, lz)
+    data = data.astype(dtype, copy=False)
     data = data[..., np.newaxis].repeat(2, axis=-1)  # Expect identical output
     with expected_warnings(['"cg" mode' + '|' + SCIPY_RANK_WARNING,
                             NUMPY_MATRIX_WARNING]):

From 6ae83c756695af4d11e2f821c848d3b4c7ef5526 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 29 Apr 2021 18:29:40 -0400
Subject: [PATCH 6/8] slic: make sure float16 gets promoted to float32

---
 skimage/segmentation/slic_superpixels.py | 11 ++++++++---
 skimage/segmentation/tests/test_slic.py  |  4 +++-
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py
index 022882f1c87..1ac99db0ab8 100644
--- a/skimage/segmentation/slic_superpixels.py
+++ b/skimage/segmentation/slic_superpixels.py
@@ -1,14 +1,16 @@
 import warnings
 from collections.abc import Iterable
+
 import numpy as np
+from numpy import random
 from scipy import ndimage as ndi
-from scipy.spatial.distance import pdist, squareform
 from scipy.cluster.vq import kmeans2
-from numpy import random
+from scipy.spatial.distance import pdist, squareform
 
-from ._slic import (_slic_cython, _enforce_label_connectivity_cython)
+from .._shared.utils import _supported_float_type
 from ..util import img_as_float, regular_grid
 from ..color import rgb2lab
+from ._slic import (_slic_cython, _enforce_label_connectivity_cython)
 
 
 def _get_mask_centroids(mask, n_centroids, multichannel):
@@ -230,6 +232,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0,
     """
 
     image = img_as_float(image)
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
     use_mask = mask is not None
     dtype = image.dtype
 
diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py
index f77cc70d798..703fdf7fdb6 100644
--- a/skimage/segmentation/tests/test_slic.py
+++ b/skimage/segmentation/tests/test_slic.py
@@ -466,7 +466,9 @@ def test_gray_3d_mask():
         assert_equal(seg[s][2:-2, 2:-2, 2:-2], c)
 
 
-@pytest.mark.parametrize("dtype", ['float32', 'float64', 'uint8', 'int'])
+@pytest.mark.parametrize(
+    "dtype", ['float16', 'float32', 'float64', 'uint8', 'int']
+)
 def test_dtype_support(dtype):
     img = np.random.rand(28, 28).astype(dtype)
 

From 0fd9c4143a079fde980cbf8b00a975ef6e7925d6 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Mon, 3 May 2021 13:39:09 -0400
Subject: [PATCH 7/8] pep8

---
 skimage/segmentation/_chan_vese.py           |  5 +++--
 skimage/segmentation/active_contour_model.py | 20 ++++++++++----------
 2 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/skimage/segmentation/_chan_vese.py b/skimage/segmentation/_chan_vese.py
index 44e8f3d2db6..f20c83ef4dd 100644
--- a/skimage/segmentation/_chan_vese.py
+++ b/skimage/segmentation/_chan_vese.py
@@ -118,7 +118,8 @@ def _cv_reset_level_set(phi):
 def _cv_checkerboard(image_size, square_size, dtype=np.float64):
     """Generates a checkerboard level set function.
 
-    According to Pascal Getreuer, such a level set function has fast convergence.
+    According to Pascal Getreuer, such a level set function has fast
+    convergence.
     """
     yv = np.arange(image_size[0], dtype=dtype).reshape(image_size[0], 1)
     xv = np.arange(image_size[1], dtype=dtype)
@@ -151,7 +152,7 @@ def _cv_small_disk(image_size):
     centerX = int((image_size[1]-1) / 2)
     res[centerY, centerX] = 0.
     radius = float(min(centerX, centerY)) / 2.0
-    return (radius - distance(res)) / (radius*3)
+    return (radius - distance(res)) / (radius * 3)
 
 
 def _cv_init_level_set(init_level_set, image_shape, dtype=np.float64):
diff --git a/skimage/segmentation/active_contour_model.py b/skimage/segmentation/active_contour_model.py
index 7f4d1e71224..65941d7761d 100644
--- a/skimage/segmentation/active_contour_model.py
+++ b/skimage/segmentation/active_contour_model.py
@@ -150,13 +150,13 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
     eye_n = np.eye(n, dtype=float)
     a = np.roll(eye_n, -1, axis=0) + \
         np.roll(eye_n, -1, axis=1) - \
-        2*eye_n  # second order derivative, central difference
+        2 * eye_n  # second order derivative, central difference
     b = np.roll(eye_n, -2, axis=0) + \
         np.roll(eye_n, -2, axis=1) - \
-        4*np.roll(eye_n, -1, axis=0) - \
-        4*np.roll(eye_n, -1, axis=1) + \
-        6*eye_n  # fourth order derivative, central difference
-    A = -alpha*a + beta*b
+        4 * np.roll(eye_n, -1, axis=0) - \
+        4 * np.roll(eye_n, -1, axis=1) + \
+        6 * eye_n  # fourth order derivative, central difference
+    A = -alpha * a + beta * b
 
     # Impose boundary conditions different from periodic:
     sfixed = False
@@ -187,7 +187,7 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
         efree = True
 
     # Only one inversion is needed for implicit spline energy minimization:
-    inv = np.linalg.inv(A + gamma*eye_n)
+    inv = np.linalg.inv(A + gamma * eye_n)
     # can use float_dtype once we have computed the inverse in double precision
     inv = inv.astype(float_dtype, copy=False)
 
@@ -213,8 +213,8 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
         yn = inv @ (gamma*y + fy)
 
         # Movements are capped to max_px_move per iteration:
-        dx = max_px_move*np.tanh(xn - x)
-        dy = max_px_move*np.tanh(yn - y)
+        dx = max_px_move * np.tanh(xn - x)
+        dy = max_px_move * np.tanh(yn - y)
         if sfixed:
             dx[0] = 0
             dy[0] = 0
@@ -231,8 +231,8 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
             xsave[j, :] = x
             ysave[j, :] = y
         else:
-            dist = np.min(np.max(np.abs(xsave - x[None, :]) +
-                                 np.abs(ysave - y[None, :]), 1))
+            dist = np.min(np.max(np.abs(xsave - x[None, :])
+                                 + np.abs(ysave - y[None, :]), 1))
             if dist < convergence:
                 break
 

From f29a1d99e2f2ecaf986decb84b1621f9ea8c4481 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <grlee77@gmail.com>
Date: Mon, 3 May 2021 21:00:02 -0400
Subject: [PATCH 8/8] Apply suggestions from code review

Co-authored-by: Riadh Fezzani <rfezzani@gmail.com>
---
 skimage/segmentation/active_contour_model.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/skimage/segmentation/active_contour_model.py b/skimage/segmentation/active_contour_model.py
index 65941d7761d..801b23c0a61 100644
--- a/skimage/segmentation/active_contour_model.py
+++ b/skimage/segmentation/active_contour_model.py
@@ -148,14 +148,14 @@ def active_contour(image, snake, alpha=0.01, beta=0.1,
 
     # Build snake shape matrix for Euler equation in double precision
     eye_n = np.eye(n, dtype=float)
-    a = np.roll(eye_n, -1, axis=0) + \
-        np.roll(eye_n, -1, axis=1) - \
-        2 * eye_n  # second order derivative, central difference
-    b = np.roll(eye_n, -2, axis=0) + \
-        np.roll(eye_n, -2, axis=1) - \
-        4 * np.roll(eye_n, -1, axis=0) - \
-        4 * np.roll(eye_n, -1, axis=1) + \
-        6 * eye_n  # fourth order derivative, central difference
+    a = (np.roll(eye_n, -1, axis=0)
+         + np.roll(eye_n, -1, axis=1)
+         - 2 * eye_n)  # second order derivative, central difference
+    b = (np.roll(eye_n, -2, axis=0)
+         + np.roll(eye_n, -2, axis=1)
+         - 4 * np.roll(eye_n, -1, axis=0)
+         - 4 * np.roll(eye_n, -1, axis=1)
+         + 6 * eye_n)  # fourth order derivative, central difference
     A = -alpha * a + beta * b
 
     # Impose boundary conditions different from periodic: