diff --git a/skimage/_shared/utils.py b/skimage/_shared/utils.py index 00f7ef3fb51..d0f3a068a9f 100644 --- a/skimage/_shared/utils.py +++ b/skimage/_shared/utils.py @@ -457,6 +457,7 @@ def _float_type(image): float32, float64, complex64, complex128 are preserved. float16 is promoted to float32. + complex256 is demoted to complex128. Other types are cast to float64. Paramters @@ -469,7 +470,11 @@ def _float_type(image): float_type : dtype Floating-point dtype for the image. """ - if image.dtype.kind in 'bui': - # case all integer types to float64 + char = image.dtype.char + if char == 'e': + return np.float32 + elif char == 'G': + return np.complex128 + elif char not in 'dfDF': return np.float64 - return np.promote_types(image.dtype, np.float32) + return image.dtype diff --git a/skimage/metrics/_structural_similarity.py b/skimage/metrics/_structural_similarity.py index 71aa363337e..f3cbc051d56 100644 --- a/skimage/metrics/_structural_similarity.py +++ b/skimage/metrics/_structural_similarity.py @@ -4,7 +4,7 @@ from ..util.dtype import dtype_range from ..util.arraycrop import crop -from .._shared.utils import warn, check_shape_equality +from .._shared.utils import _float_type, warn, check_shape_equality __all__ = ['structural_similarity'] @@ -87,6 +87,7 @@ def structural_similarity(im1, im2, """ check_shape_equality(im1, im2) + float_type = _float_type(im1) if multichannel: # loop over channels @@ -98,11 +99,11 @@ def structural_similarity(im1, im2, full=full) args.update(kwargs) nch = im1.shape[-1] - 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) for ch in range(nch): ch_result = structural_similarity(im1[..., ch], im2[..., ch], **args) @@ -173,8 +174,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 @@ -210,8 +211,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. diff --git a/skimage/metrics/simple_metrics.py b/skimage/metrics/simple_metrics.py index e4474b5d72b..968ce849770 100644 --- a/skimage/metrics/simple_metrics.py +++ b/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 _float_type, warn, check_shape_equality __all__ = ['mean_squared_error', 'normalized_root_mse', diff --git a/skimage/metrics/tests/test_simple_metrics.py b/skimage/metrics/tests/test_simple_metrics.py index 1d90a48f596..df8e9859c27 100644 --- a/skimage/metrics/tests/test_simple_metrics.py +++ b/skimage/metrics/tests/test_simple_metrics.py @@ -1,5 +1,6 @@ from skimage._shared._warnings import expected_warnings from skimage._shared.testing import assert_equal, assert_almost_equal +from skimage._shared.utils import _float_type from skimage._shared import testing import numpy as np @@ -38,10 +39,13 @@ def test_PSNR_vs_IPOL(): assert_almost_equal(p, p_IPOL, decimal=4) -def test_PSNR_float(): +@testing.parametrize('dtype', [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) + 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 assert_almost_equal(p_uint8, p_float64, decimal=5) # mixed precision inputs @@ -62,11 +66,13 @@ def test_PSNR_errors(): 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)) +@testing.parametrize('dtype', [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'), @@ -107,14 +113,13 @@ def test_nmi_different_sizes(): assert normalized_mutual_information(cam[:, :400], cam[:400, :]) > 1 -def test_nmi_random(): +@testing.parametrize('dtype', [np.float32, np.float64]) +def test_nmi_random(dtype): 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, - ) + 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(): diff --git a/skimage/metrics/tests/test_structural_similarity.py b/skimage/metrics/tests/test_structural_similarity.py index 692526e8932..d6ea3a910c4 100644 --- a/skimage/metrics/tests/test_structural_similarity.py +++ b/skimage/metrics/tests/test_structural_similarity.py @@ -8,6 +8,7 @@ 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 _float_type np.random.seed(5) cam = data.camera() @@ -53,7 +54,8 @@ 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): +@testing.parametrize('dtype', [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,15 +79,19 @@ def test_structural_similarity_grad(seed): mssim, grad, s = structural_similarity( X, Y, data_range=255, gradient=True, full=True) + s.dtype == dtype + grad.dtype == dtype assert np.all(grad < 0.05) -def test_structural_similarity_dtype(): +@testing.parametrize('dtype', [np.float32, np.float64]) +def test_structural_similarity_dtype(dtype): N = 30 - X = np.random.rand(N, N) - Y = np.random.rand(N, N) + X = np.random.rand(N, N).astype(dtype, copy=False) + Y = np.random.rand(N, N).astype(dtype, copy=False) S1 = structural_similarity(X, Y) + assert S1.dtype == np.float64 X = (X * 255).astype(np.uint8) Y = (X * 255).astype(np.uint8) @@ -128,15 +134,17 @@ def test_structural_similarity_multichannel(): structural_similarity(Xc, Yc, win_size=7, multichannel=False) -def test_structural_similarity_nD(): +@testing.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