Promote small integer types to single rather than double precision

python/cucim/src/cucim/skimage/_shared/tests/test_utils.py

@@ -157,7 +157,8 @@ def test_validate_interpolation_order(dtype, order):
 )
 def test_supported_float_dtype_real(dtype):
     float_dtype = _supported_float_type(dtype)
-    if dtype in [np.float16, np.float32]:
+    if dtype in [np.float16, np.float32, np.int8, np.uint8, np.int16,
+                 np.uint16, bool]:
         assert float_dtype == np.float32
     else:
         assert float_dtype == np.float64
@@ -188,7 +189,9 @@ def test_supported_float_dtype_input_kinds(dtype):
     'dtypes, expected',
     [
         ((np.float16, np.float64), np.float64),
-        ([np.float32, np.uint16, np.int8], np.float64),
+        ([np.float32, np.uint16, np.int8], np.float32),
+        ([np.float32, bool], np.float32),
+        ([np.float32, np.uint32, np.int16], np.float64),
         ({np.float32, np.float16}, np.float32),
     ]
 )

python/cucim/src/cucim/skimage/_shared/utils.py

@@ -654,9 +654,10 @@ def convert_to_float(image, preserve_range):
         # Convert image to double only if it is not single or double
         # precision float
         if image.dtype.char not in 'df':
-            image = image.astype(float)
+            image = image.astype(_supported_float_type(image.dtype))
     else:
         from ..util.dtype import img_as_float
+
         image = img_as_float(image)
     return image
 
@@ -730,14 +731,21 @@ def _fix_ndimage_mode(mode):
 
 new_float_type = {
     # preserved types
-    cp.float32().dtype.char: cp.float32,
-    cp.float64().dtype.char: cp.float64,
-    cp.complex64().dtype.char: cp.complex64,
-    cp.complex128().dtype.char: cp.complex128,
-    # altered types
-    cp.float16().dtype.char: cp.float32,
-    'g': cp.float64,      # cp.float128 ; doesn't exist on windows
-    'G': cp.complex128,   # cp.complex256 ; doesn't exist on windows
+    'f': cp.float32,     # float32
+    'd': cp.float64,     # float64
+    'F': cp.complex64,   # complex64
+    'D': cp.complex128,  # complex128
+    # promoted float types
+    'e': cp.float32,     # float16
+    # truncated float types
+    'g': cp.float64,     # float128 (doesn't exist on windows)
+    'G': cp.complex128,  # complex256 (doesn't exist on windows)
+    # integer types that can be exactly represented in float32
+    'b': cp.float32,     # int8
+    'B': cp.float32,     # uint8
+    'h': cp.float32,     # int16
+    'H': cp.float32,     # uint16
+    '?': cp.float32,     # bool
 }
 
 

python/cucim/src/cucim/skimage/color/tests/test_colorconv.py

@@ -14,7 +14,8 @@
 import cupy as cp
 import numpy as np
 import pytest
-from cupy.testing import assert_array_almost_equal, assert_array_equal
+from cupy.testing import (assert_allclose, assert_array_almost_equal,
+                          assert_array_equal)
 from numpy.testing import assert_equal
 from skimage import data
 
@@ -236,7 +237,7 @@ def test_xyz_rgb_roundtrip(self, channel_axis):
         round_trip = xyz2rgb(rgb2xyz(img_rgb, channel_axis=channel_axis),
                              channel_axis=channel_axis)
 
-        assert_array_almost_equal(round_trip, img_rgb)
+        assert_allclose(round_trip, img_rgb, rtol=1e-5, atol=1e-5)
 
     # RGB<->HED roundtrip with ubyte image
     def test_hed_rgb_roundtrip(self):
@@ -497,12 +498,15 @@ def test_rgb2lab_brucelindbloom(self):
     def test_lab_rgb_roundtrip(self, channel_axis):
         img_rgb = img_as_float(self.img_rgb)
         img_rgb = cp.moveaxis(img_rgb, source=-1, destination=channel_axis)
-        assert_array_almost_equal(
+
+        assert_allclose(
             lab2rgb(
                 rgb2lab(img_rgb, channel_axis=channel_axis),
                 channel_axis=channel_axis
             ),
             img_rgb,
+            rtol=1e-5,
+            atol=1e-5,
         )
 
     def test_rgb2lab_dtype(self):
@@ -633,12 +637,14 @@ def test_luv2rgb_dtype(self):
     def test_luv_rgb_roundtrip(self, channel_axis):
         img_rgb = img_as_float(self.img_rgb)
         img_rgb = cp.moveaxis(img_rgb, source=-1, destination=channel_axis)
-        assert_array_almost_equal(
+        assert_allclose(
             luv2rgb(
                 rgb2luv(img_rgb, channel_axis=channel_axis),
                 channel_axis=channel_axis
             ),
             img_rgb,
+            rtol=1e-4,
+            atol=1e-4,
         )
 
     def test_lab_rgb_outlier(self):
@@ -674,15 +680,15 @@ def test_lab_lch_roundtrip(self, channel_axis):
             lab2lch(lab, channel_axis=channel_axis),
             channel_axis=channel_axis,
         )
-        assert_array_almost_equal(lab2, lab)
+        assert_allclose(lab2, lab, rtol=1e-4, atol=1e-4)
 
     def test_rgb_lch_roundtrip(self):
         rgb = img_as_float(self.img_rgb)
         lab = rgb2lab(rgb)
         lch = lab2lch(lab)
         lab2 = lch2lab(lch)
         rgb2 = lab2rgb(lab2)
-        assert_array_almost_equal(rgb, rgb2)
+        assert_allclose(rgb, rgb2, rtol=1e-4, atol=1e-4)
 
     def test_lab_lch_0d(self):
         lab0 = self._get_lab0()
@@ -736,26 +742,41 @@ def test_yuv(self):
     def test_yuv_roundtrip(self, channel_axis):
         img_rgb = img_as_float(self.img_rgb)[::16, ::16]
         img_rgb = cp.moveaxis(img_rgb, source=-1, destination=channel_axis)
-        assert_array_almost_equal(
+        assert_allclose(
             yuv2rgb(rgb2yuv(img_rgb, channel_axis=channel_axis),
                     channel_axis=channel_axis),
-            img_rgb)
-        assert_array_almost_equal(
+            img_rgb,
+            rtol=1e-5,
+            atol=1e-5,
+        )
+        assert_allclose(
             yiq2rgb(rgb2yiq(img_rgb, channel_axis=channel_axis),
                     channel_axis=channel_axis),
-            img_rgb)
-        assert_array_almost_equal(
+            img_rgb,
+            rtol=1e-5,
+            atol=1e-5,
+        )
+        assert_allclose(
             ypbpr2rgb(rgb2ypbpr(img_rgb, channel_axis=channel_axis),
                       channel_axis=channel_axis),
-            img_rgb)
-        assert_array_almost_equal(
+            img_rgb,
+            rtol=1e-5,
+            atol=1e-5,
+        )
+        assert_allclose(
             ycbcr2rgb(rgb2ycbcr(img_rgb, channel_axis=channel_axis),
                       channel_axis=channel_axis),
-            img_rgb)
-        assert_array_almost_equal(
+            img_rgb,
+            rtol=1e-5,
+            atol=1e-5,
+        )
+        assert_allclose(
             ydbdr2rgb(rgb2ydbdr(img_rgb, channel_axis=channel_axis),
                       channel_axis=channel_axis),
-            img_rgb)
+            img_rgb,
+            rtol=1e-5,
+            atol=1e-5,
+        )
 
     def test_rgb2yuv_dtype(self):
         img = self.colbars_array.astype('float64')

python/cucim/src/cucim/skimage/exposure/exposure.py

@@ -440,7 +440,7 @@ def intensity_range(image, range_values="image", clip_negative=False):
     return i_min, i_max
 
 
-def _output_dtype(dtype_or_range):
+def _output_dtype(dtype_or_range, image_dtype):
     """Determine the output dtype for rescale_intensity.
 
     The dtype is determined according to the following rules:
@@ -450,13 +450,16 @@ def _output_dtype(dtype_or_range):
       in which case the data type that can contain it will be used
       (e.g. uint16 in this case).
     - if ``dtype_or_range`` is a pair of values, the output data type will be
-      float.
+      ``_supported_float_type(image_dtype)``. This preserves float32 output for
+      float32 inputs.
 
     Parameters
     ----------
     dtype_or_range : type, string, or 2-tuple of int/float
         The desired range for the output, expressed as either a NumPy dtype or
         as a (min, max) pair of numbers.
+    image_dtype : np.dtype
+        The input image dtype.
 
     Returns
     -------
@@ -465,7 +468,7 @@ def _output_dtype(dtype_or_range):
     """
     if type(dtype_or_range) in [list, tuple, np.ndarray]:
         # pair of values: always return float.
-        return float
+        return utils._supported_float_type(image_dtype)
     if type(dtype_or_range) == type:
         # already a type: return it
         return dtype_or_range
@@ -577,9 +580,9 @@ def rescale_intensity(image, in_range="image", out_range="dtype"):
     array([127, 127, 127], dtype=int32)
     """
     if out_range in ['dtype', 'image']:
-        out_dtype = _output_dtype(image.dtype.type)
+        out_dtype = _output_dtype(image.dtype.type, image.dtype)
     else:
-        out_dtype = _output_dtype(out_range)
+        out_dtype = _output_dtype(out_range, image.dtype)
 
     imin, imax = map(float, intensity_range(image, in_range))
     omin, omax = map(float, intensity_range(image, out_range,

python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py

@@ -387,7 +387,7 @@ def test_rescale_float_output():
     image = cp.array([-128, 0, 127], dtype=cp.int8)
     output_image = exposure.rescale_intensity(image, out_range=(0, 255))
     cp.testing.assert_array_equal(output_image, [0, 128, 255])
-    assert output_image.dtype == float
+    assert output_image.dtype == _supported_float_type(image.dtype)
 
 
 def test_rescale_raises_on_incorrect_out_range():

python/cucim/src/cucim/skimage/feature/corner.py

@@ -140,8 +140,13 @@ def structure_tensor(image, sigma=1, mode="constant", cval=0, order=None):
     if order == "xy":
         derivatives = reversed(derivatives)
 
+    # Autodetection as done internally to Gaussian, but set it here to silence
+    # a warning.
+    channel_axis = -1 if (image.ndim == 3 and image.shape[-1] == 3) else None
+
     # structure tensor
-    A_elems = [gaussian(der0 * der1, sigma, mode=mode, cval=cval)
+    A_elems = [gaussian(der0 * der1, sigma, mode=mode, cval=cval,
+                        channel_axis=channel_axis)
                for der0, der1 in combinations_with_replacement(derivatives, 2)]
 
     return A_elems
@@ -205,7 +210,12 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0, order='rc'):
     float_dtype = _supported_float_type(image.dtype)
     image = image.astype(float_dtype, copy=False)
 
-    gaussian_filtered = gaussian(image, sigma=sigma, mode=mode, cval=cval)
+    # Autodetection as done internally to Gaussian, but set it here to silence
+    # a warning.
+    channel_axis = -1 if (image.ndim == 3 and image.shape[-1] == 3) else None
+
+    gaussian_filtered = gaussian(image, sigma=sigma, mode=mode, cval=cval,
+                                 channel_axis=channel_axis)
 
     gradients = cp.gradient(gaussian_filtered)
     axes = range(image.ndim)

python/cucim/src/cucim/skimage/feature/template.py

@@ -134,7 +134,9 @@ def match_template(image, template, pad_input=False, mode='constant',
     image_shape = image.shape
 
     float_dtype = _supported_float_type(image.dtype)
-    image = image.astype(float_dtype, copy=False)
+
+    # Note: keep image in float64 for accuracy of cumsum operations, etc.
+    image = image.astype(cp.float64, copy=False)
     template = template.astype(float_dtype, copy=False)
 
     pad_width = tuple((width, width) for width in template.shape)
@@ -153,11 +155,12 @@ def match_template(image, template, pad_input=False, mode='constant',
         image_window_sum = _window_sum_3d(image, template.shape)
         image_window_sum2 = _window_sum_3d(image * image, template.shape)
 
-    template_mean = template.mean()
+    # perform mean and sum in float64 for accuracy
+    template_mean = template.mean(dtype=cp.float64)
     template_volume = _misc.prod(template.shape)
     template_ssd = template - template_mean
     template_ssd *= template_ssd
-    template_ssd = cp.sum(template_ssd)
+    template_ssd = cp.sum(template_ssd, dtype=cp.float64)
 
     if image.ndim == 2:
         xcorr = signal.fftconvolve(image, template[::-1, ::-1],

python/cucim/src/cucim/skimage/feature/tests/test_template.py

@@ -26,7 +26,7 @@ def test_template(dtype):
     target = cp.asarray(target)
 
     result = match_template(image, target)
-    assert result.dtype == dtype
+    assert result.dtype == result.dtype
     delta = 5
 
     positions = peak_local_max(result, min_distance=delta)
@@ -200,6 +200,5 @@ def test_bounding_values():
     template = cp.zeros((3, 3))
     template[1, 1] = 1
     result = match_template(image, template)
-    print(result.max())
     assert result.max() < 1 + 1e-7
     assert result.min() > -1 - 1e-7

python/cucim/src/cucim/skimage/filters/_fft_based.py

@@ -123,6 +123,7 @@ def butterworth(
                  else np.delete(image.shape, channel_axis))
     is_real = cp.isrealobj(image)
     float_dtype = _supported_float_type(image.dtype, allow_complex=True)
+    image = image.astype(float_dtype, copy=False)
     wfilt = _get_ND_butterworth_filter(
         fft_shape, cutoff_frequency_ratio, order, high_pass, is_real,
         float_dtype

python/cucim/src/cucim/skimage/filters/ridges.py

@@ -234,13 +234,13 @@ def meijering(image, sigmas=range(1, 10, 2), alpha=None,
     if alpha is None:
         alpha = 1.0 / ndim
 
-    float_dtype = _supported_float_type(image.dtype)
-    image = image.astype(float_dtype, copy=False)
-
     # Invert image to detect dark ridges on bright background
     if black_ridges:
         image = invert(image)
 
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
     # Generate empty (n+1)D arrays for storing auxiliary images filtered at
     # different (sigma) scales
     filtered_array = cp.empty(sigmas.shape + image.shape, dtype=float_dtype)

python/cucim/src/cucim/skimage/filters/tests/test_fft_based.py

@@ -108,7 +108,7 @@ def test_butterworth_4D_channel(chan, dtype):
 
 
 def test_butterworth_correctness_bw():
-    small = cp.array(coins()[180:190, 260:270])
+    small = cp.array(coins()[180:190, 260:270], dtype=float)
     filtered = butterworth(small,
                            cutoff_frequency_ratio=0.2)
     correct = cp.array(
@@ -129,7 +129,7 @@ def test_butterworth_correctness_bw():
 
 
 def test_butterworth_correctness_rgb():
-    small = cp.array(astronaut()[135:145, 205:215])
+    small = cp.array(astronaut()[135:145, 205:215], dtype=float)
     filtered = butterworth(small,
                            cutoff_frequency_ratio=0.3,
                            high_pass=True,

python/cucim/src/cucim/skimage/filters/tests/test_ridges.py

@@ -4,7 +4,7 @@
 from cupy.testing import assert_allclose, assert_array_equal, assert_array_less
 from skimage.data import camera, retina
 
-from cucim.skimage import img_as_float
+from cucim.skimage import img_as_float, img_as_float64
 from cucim.skimage._shared.utils import _supported_float_type
 from cucim.skimage.color import rgb2gray
 from cucim.skimage.filters import frangi, hessian, meijering, sato
@@ -178,9 +178,12 @@ def test_3d_linearity():
                     atol=1e-3)
 
 
-def test_2d_cropped_camera_image():
-
+@pytest.mark.parametrize('dtype', ['float64', 'uint8'])
+def test_2d_cropped_camera_image(dtype):
     a_black = crop(cp.array(camera()), ((200, 212), (100, 312)))
+    assert a_black.dtype == cp.uint8
+    if dtype == 'float64':
+        a_black = img_as_float64(a_black)
     a_white = invert(a_black)
 
     zeros = cp.zeros((100, 100))
@@ -208,9 +211,13 @@ def test_ridge_output_dtype(func, dtype):
     assert func(img).dtype == _supported_float_type(img.dtype)
 
 
-def test_3d_cropped_camera_image():
+@pytest.mark.parametrize('dtype', ['float64', 'uint8'])
+def test_3d_cropped_camera_image(dtype):
 
     a_black = crop(cp.asarray(camera()), ((200, 212), (100, 312)))
+    assert a_black.dtype == cp.uint8
+    if dtype == 'float64':
+        a_black = img_as_float64(a_black)
     a_black = cp.dstack([a_black, a_black, a_black])
     a_white = invert(a_black)