single precision support in skimage.metrics module

skimage/metrics/_structural_similarity.py

@@ -1,13 +1,12 @@
 import functools
-from warnings import warn
 
 import numpy as np
 from scipy.ndimage import uniform_filter, gaussian_filter
 
+from .._shared import utils
+from .._shared.utils import _supported_float_type, check_shape_equality, warn
 from ..util.dtype import dtype_range
 from ..util.arraycrop import crop
-from .._shared import utils
-from .._shared.utils import warn, check_shape_equality
 
 
 __all__ = ['structural_similarity']
@@ -100,6 +99,7 @@ def structural_similarity(im1, im2,
 
     """
     check_shape_equality(im1, im2)
+    float_type = _supported_float_type(im1.dtype)
 
     if channel_axis is not None:
         # loop over channels
@@ -111,11 +111,12 @@ def structural_similarity(im1, im2,
                     full=full)
         args.update(kwargs)
         nch = im1.shape[channel_axis]
-        mssim = np.empty(nch)
+        mssim = np.empty(nch, dtype=float_type)
+
         if gradient:
-            G = np.empty(im1.shape)
+            G = np.empty(im1.shape, dtype=float_type)
         if full:
-            S = np.empty(im1.shape)
+            S = np.empty(im1.shape, dtype=float_type)
         channel_axis = channel_axis % im1.ndim
         _at = functools.partial(utils.slice_at_axis, axis=channel_axis)
         for ch in range(nch):
@@ -189,8 +190,8 @@ def structural_similarity(im1, im2,
         filter_args = {'size': win_size}
 
     # ndimage filters need floating point data
-    im1 = im1.astype(np.float64)
-    im2 = im2.astype(np.float64)
+    im1 = im1.astype(float_type, copy=False)
+    im2 = im2.astype(float_type, copy=False)
 
     NP = win_size ** ndim
 
@@ -226,8 +227,8 @@ def structural_similarity(im1, im2,
     # to avoid edge effects will ignore filter radius strip around edges
     pad = (win_size - 1) // 2
 
-    # compute (weighted) mean of ssim
-    mssim = crop(S, pad).mean()
+    # compute (weighted) mean of ssim. Use float64 for accuracy.
+    mssim = crop(S, pad).mean(dtype=np.float64)
 
     if gradient:
         # The following is Eqs. 7-8 of Avanaki 2009.

skimage/metrics/simple_metrics.py

@@ -2,7 +2,7 @@
 from scipy.stats import entropy
 
 from ..util.dtype import dtype_range
-from .._shared.utils import warn, check_shape_equality
+from .._shared.utils import _supported_float_type, check_shape_equality, warn
 
 __all__ = ['mean_squared_error',
            'normalized_root_mse',
@@ -15,7 +15,7 @@ def _as_floats(image0, image1):
     """
     Promote im1, im2 to nearest appropriate floating point precision.
     """
-    float_type = np.result_type(image0.dtype, image1.dtype, np.float32)
+    float_type = _supported_float_type([image0.dtype, image1.dtype])
     image0 = np.asarray(image0, dtype=float_type)
     image1 = np.asarray(image1, dtype=float_type)
     return image0, image1

skimage/metrics/tests/test_set_metrics.py

@@ -1,12 +1,10 @@
-from __future__ import print_function, division
+import itertools
 
 import numpy as np
+import pytest
 from numpy.testing import assert_almost_equal, assert_array_equal
 from scipy.spatial import distance
-import itertools
-import pytest
 
-from skimage._shared.testing import parametrize
 from skimage.metrics import hausdorff_distance, hausdorff_pair
 
 
@@ -37,8 +35,10 @@ def test_hausdorff_simple():
                                                             points_b))
 
 
-@parametrize("points_a, points_b",
-             itertools.product([(0, 0), (3, 0), (1, 4), (4, 1)], repeat=2))
+@pytest.mark.parametrize(
+    "points_a, points_b",
+    itertools.product([(0, 0), (3, 0), (1, 4), (4, 1)], repeat=2)
+)
 def test_hausdorff_region_single(points_a, points_b):
     shape = (5, 5)
     coords_a = np.zeros(shape, dtype=bool)
@@ -53,9 +53,11 @@ def test_hausdorff_region_single(points_a, points_b):
                                                             points_b))
 
 
-@parametrize("points_a, points_b",
-             itertools.product([(5, 4), (4, 5), (3, 4), (4, 3)],
-                               [(6, 4), (2, 6), (2, 4), (4, 0)]))
+@pytest.mark.parametrize(
+    "points_a, points_b",
+    itertools.product([(5, 4), (4, 5), (3, 4), (4, 3)],
+                      [(6, 4), (2, 6), (2, 4), (4, 0)])
+)
 def test_hausdorff_region_different_points(points_a, points_b):
     shape = (7, 7)
     coords_a = np.zeros(shape, dtype=bool)
@@ -118,9 +120,11 @@ def test_gallery():
            np.equal(hd_points, ((30, 40), (30, 50))).all()
 
 
-@parametrize("points_a, points_b",
-             itertools.product([(0, 0, 1), (0, 1, 0), (1, 0, 0)],
-                               [(0, 0, 2), (0, 2, 0), (2, 0, 0)]))
+@pytest.mark.parametrize(
+    "points_a, points_b",
+    itertools.product([(0, 0, 1), (0, 1, 0), (1, 0, 0)],
+                      [(0, 0, 2), (0, 2, 0), (2, 0, 0)])
+)
 def test_3d_hausdorff_region(points_a, points_b):
     shape = (3, 3, 3)
     coords_a = np.zeros(shape, dtype=bool)

skimage/metrics/tests/test_simple_metrics.py

@@ -1,9 +1,10 @@
-from skimage._shared._warnings import expected_warnings
-from skimage._shared.testing import assert_equal, assert_almost_equal
-from skimage._shared import testing
 import numpy as np
+import pytest
+from numpy.testing import assert_equal, assert_almost_equal
 
 from skimage import data
+from skimage._shared._warnings import expected_warnings
+from skimage._shared.utils import _supported_float_type
 from skimage.metrics import (peak_signal_noise_ratio, normalized_root_mse,
                              mean_squared_error, normalized_mutual_information)
 
@@ -38,35 +39,48 @@ def test_PSNR_vs_IPOL():
     assert_almost_equal(p, p_IPOL, decimal=4)
 
 
-def test_PSNR_float():
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_PSNR_float(dtype):
     p_uint8 = peak_signal_noise_ratio(cam, cam_noisy)
-    p_float64 = peak_signal_noise_ratio(cam / 255., cam_noisy / 255.,
-                                        data_range=1)
-    assert_almost_equal(p_uint8, p_float64, decimal=5)
+    camf = (cam / 255.).astype(dtype, copy=False)
+    camf_noisy = (cam_noisy / 255.).astype(dtype, copy=False)
+    p_float64 = peak_signal_noise_ratio(camf, camf_noisy, data_range=1)
+    assert p_float64.dtype == np.float64
+    decimal = 3 if dtype == np.float16 else 5
+    assert_almost_equal(p_uint8, p_float64, decimal=decimal)
 
     # mixed precision inputs
     p_mixed = peak_signal_noise_ratio(cam / 255., np.float32(cam_noisy / 255.),
                                       data_range=1)
-    assert_almost_equal(p_mixed, p_float64, decimal=5)
+
+    assert_almost_equal(p_mixed, p_float64, decimal=decimal)
+
+    # mismatched dtype results in a warning if data_range is unspecified
+    with expected_warnings(['Inputs have mismatched dtype']):
+        p_mixed = peak_signal_noise_ratio(cam / 255.,
+                                          np.float32(cam_noisy / 255.))
+    assert_almost_equal(p_mixed, p_float64, decimal=decimal)
 
     # mismatched dtype results in a warning if data_range is unspecified
     with expected_warnings(['Inputs have mismatched dtype']):
         p_mixed = peak_signal_noise_ratio(cam / 255.,
                                           np.float32(cam_noisy / 255.))
-    assert_almost_equal(p_mixed, p_float64, decimal=5)
+    assert_almost_equal(p_mixed, p_float64, decimal=decimal)
 
 
 def test_PSNR_errors():
     # shape mismatch
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         peak_signal_noise_ratio(cam, cam[:-1, :])
 
 
-def test_NRMSE():
-    x = np.ones(4)
-    y = np.asarray([0., 2., 2., 2.])
-    assert_equal(normalized_root_mse(y, x, normalization='mean'),
-                 1 / np.mean(y))
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_NRMSE(dtype):
+    x = np.ones(4, dtype=dtype)
+    y = np.asarray([0., 2., 2., 2.], dtype=dtype)
+    nrmse = normalized_root_mse(y, x, normalization='mean')
+    assert nrmse.dtype == np.float64
+    assert_equal(nrmse, 1 / np.mean(y))
     assert_equal(normalized_root_mse(y, x, normalization='euclidean'),
                  1 / np.sqrt(3))
     assert_equal(normalized_root_mse(y, x, normalization='min-max'),
@@ -90,10 +104,10 @@ def test_NRMSE_no_int_overflow():
 def test_NRMSE_errors():
     x = np.ones(4)
     # shape mismatch
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         normalized_root_mse(x[:-1], x)
     # invalid normalization name
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         normalized_root_mse(x, x, normalization='foo')
 
 
@@ -107,14 +121,14 @@ def test_nmi_different_sizes():
     assert normalized_mutual_information(cam[:, :400], cam[:400, :]) > 1
 
 
-def test_nmi_random():
-    random1 = np.random.random((100, 100))
-    random2 = np.random.random((100, 100))
-    assert_almost_equal(
-        normalized_mutual_information(random1, random2, bins=10),
-        1,
-        decimal=2,
-    )
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_nmi_random(dtype):
+    rng = np.random.default_rng()
+    random1 = rng.random((100, 100)).astype(dtype)
+    random2 = rng.random((100, 100)).astype(dtype)
+    nmi = normalized_mutual_information(random1, random2, bins=10)
+    assert nmi.dtype == np.float64
+    assert_almost_equal(nmi, 1, decimal=2)
 
 
 def test_nmi_random_3d():

skimage/metrics/tests/test_structural_similarity.py

@@ -1,13 +1,14 @@
 import os
+
 import numpy as np
+import pytest
 
 from skimage import data
-from skimage.metrics import structural_similarity
-
-from skimage._shared import testing
 from skimage._shared._warnings import expected_warnings
 from skimage._shared.testing import (assert_equal, assert_almost_equal,
                                      assert_array_almost_equal, fetch)
+from skimage._shared.utils import _supported_float_type
+from skimage.metrics import structural_similarity
 
 np.random.seed(5)
 cam = data.camera()
@@ -52,8 +53,9 @@ def test_structural_similarity_image():
 
 # Because we are forcing a random seed state, it is probably good to test
 # against a few seeds in case on seed gives a particularly bad example
-@testing.parametrize('seed', [1, 2, 3, 5, 8, 13])
-def test_structural_similarity_grad(seed):
+@pytest.mark.parametrize('seed', [1, 2, 3, 5, 8, 13])
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
+def test_structural_similarity_grad(seed, dtype):
     N = 30
     # NOTE: This test is known to randomly fail on some systems (Mac OS X 10.6)
     #       And when testing tests in parallel. Therefore, we choose a few
@@ -64,8 +66,8 @@ def test_structural_similarity_grad(seed):
     # X = np.random.rand(N, N) * 255
     # Y = np.random.rand(N, N) * 255
     rnd = np.random.RandomState(seed)
-    X = rnd.rand(N, N) * 255
-    Y = rnd.rand(N, N) * 255
+    X = rnd.rand(N, N).astype(dtype, copy=False) * 255
+    Y = rnd.rand(N, N).astype(dtype, copy=False) * 255
 
     f = structural_similarity(X, Y, data_range=255)
     g = structural_similarity(X, Y, data_range=255, gradient=True)
@@ -77,26 +79,32 @@ def test_structural_similarity_grad(seed):
 
     mssim, grad, s = structural_similarity(
         X, Y, data_range=255, gradient=True, full=True)
+    assert s.dtype == _supported_float_type(dtype)
+    assert grad.dtype == _supported_float_type(dtype)
     assert np.all(grad < 0.05)
 
 
-def test_structural_similarity_dtype():
+@pytest.mark.parametrize(
+    'dtype', [np.uint8, np.int32, np.float16, np.float32, np.float64]
+)
+def test_structural_similarity_dtype(dtype):
     N = 30
     X = np.random.rand(N, N)
     Y = np.random.rand(N, N)
+    if np.dtype(dtype).kind in 'iub':
+        X = (X * 255).astype(np.uint8)
+        Y = (X * 255).astype(np.uint8)
+    else:
+        X = X.astype(dtype, copy=False)
+        Y = Y.astype(dtype, copy=False)
 
     S1 = structural_similarity(X, Y)
-
-    X = (X * 255).astype(np.uint8)
-    Y = (X * 255).astype(np.uint8)
-
-    S2 = structural_similarity(X, Y)
+    assert S1.dtype == np.float64
 
     assert S1 < 0.1
-    assert S2 < 0.1
 
 
-@testing.parametrize('channel_axis', [0, 1, 2, -1])
+@pytest.mark.parametrize('channel_axis', [0, 1, 2, -1])
 def test_structural_similarity_multichannel(channel_axis):
     N = 100
     X = (np.random.rand(N, N) * 255).astype(np.uint8)
@@ -132,7 +140,7 @@ def test_structural_similarity_multichannel(channel_axis):
     assert_equal(S3.shape, Xc.shape)
 
     # fail if win_size exceeds any non-channel dimension
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         structural_similarity(Xc, Yc, win_size=7, channel_axis=None)
 
 
@@ -151,15 +159,17 @@ def test_structural_similarity_multichannel_deprecated():
     assert_almost_equal(S1, S2)
 
 
-def test_structural_similarity_nD():
+@pytest.mark.parametrize('dtype', [np.uint8, np.float32, np.float64])
+def test_structural_similarity_nD(dtype):
     # test 1D through 4D on small random arrays
     N = 10
     for ndim in range(1, 5):
         xsize = [N, ] * 5
-        X = (np.random.rand(*xsize) * 255).astype(np.uint8)
-        Y = (np.random.rand(*xsize) * 255).astype(np.uint8)
+        X = (np.random.rand(*xsize) * 255).astype(dtype)
+        Y = (np.random.rand(*xsize) * 255).astype(dtype)
 
         mssim = structural_similarity(X, Y, win_size=3)
+        assert mssim.dtype == np.float64
         assert mssim < 0.05
 
 
@@ -227,15 +237,15 @@ def test_invalid_input():
     # size mismatch
     X = np.zeros((9, 9), dtype=np.double)
     Y = np.zeros((8, 8), dtype=np.double)
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         structural_similarity(X, Y)
     # win_size exceeds image extent
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         structural_similarity(X, X, win_size=X.shape[0] + 1)
     # some kwarg inputs must be non-negative
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         structural_similarity(X, X, K1=-0.1)
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         structural_similarity(X, X, K2=-0.1)
-    with testing.raises(ValueError):
+    with pytest.raises(ValueError):
         structural_similarity(X, X, sigma=-1.0)