single precision support in skimage.segmentation

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,18 @@ 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.
+    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 +139,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 +152,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 +169,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 +302,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)

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.")

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,32 +44,40 @@ 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.
     # Reference implementation doesn't do that, though, and it only has
     # 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:

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,21 +140,23 @@ 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
-    A = -alpha*a + beta*b
+    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:
     sfixed = False
@@ -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
 

skimage/segmentation/boundaries.py

@@ -1,6 +1,8 @@
 
 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 +214,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':

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