diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 9e8df953..be5f2bd0 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -2,6 +2,7 @@ from ._find_contours import find_contours from ._marching_cubes import (marching_cubes, mesh_surface_area, correct_mesh_orientation) from ._regionprops import regionprops, perimeter +from .simple_metrics import mean_squared_error, normalized_root_mse, psnr from ._structural_similarity import structural_similarity from ._polygon import approximate_polygon, subdivide_polygon from ._pnpoly import points_in_poly, grid_points_in_poly @@ -34,4 +35,7 @@ __all__ = ['find_contours', 'profile_line', 'label', 'points_in_poly', - 'grid_points_in_poly'] + 'grid_points_in_poly', + 'mean_squared_error', + 'normalized_root_mse', + 'psnr'] diff --git a/skimage/measure/simple_metrics.py b/skimage/measure/simple_metrics.py new file mode 100644 index 00000000..6083010b --- /dev/null +++ b/skimage/measure/simple_metrics.py @@ -0,0 +1,132 @@ +from __future__ import division + +import numpy as np +from ..util.dtype import dtype_range + +__all__ = ['mean_squared_error', 'normalized_root_mse', 'psnr'] + + +def _assert_compatible(im1, im2): + """Raise an error if the shape and dtype do not match.""" + if not im1.dtype == im2.dtype: + raise ValueError('Input images must have the same dtype.') + if not im1.shape == im2.shape: + raise ValueError('Input images must have the same dimensions.') + return + + +def _as_floats(im1, im2): + """Promote im1, im2 to nearest appropriate floating point precision.""" + float_type = np.result_type(im1.dtype, im2.dtype, np.float32) + if im1.dtype != float_type: + im1 = im1.astype(float_type) + if im2.dtype != float_type: + im2 = im2.astype(float_type) + return im1, im2 + + +def mean_squared_error(im1, im2): + """Compute the mean-squared error between two images. + + Parameters + ---------- + im1, im2 : ndarray + Image. Any dimensionality. + + Returns + ------- + mse : float + The mean-squared error (MSE) metric. + + """ + _assert_compatible(im1, im2) + im1, im2 = _as_floats(im1, im2) + return np.mean(np.square(im1 - im2), dtype=np.float64) + + +def normalized_root_mse(im_true, im_test, norm_type='Euclidean'): + """Compute the normalized root mean-squared error (NRMSE) between two + images. + + Parameters + ---------- + im_true : ndarray + Ground-truth image. + im_test : ndarray + Test image. + norm_type : {'Euclidean', 'min-max', 'mean'} + Controls the normalization method to use in the denominator of the + NRMSE. There is no standard method of normalization across the + literature [1]_. The methods available here are as follows: + + - 'Euclidean' : normalize by the Euclidean norm of ``im_true``. + - 'min-max' : normalize by the intensity range of ``im_true``. + - 'mean' : normalize by the mean of ``im_true``. + + Returns + ------- + nrmse : float + The NRMSE metric. + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Root-mean-square_deviation + + """ + _assert_compatible(im_true, im_test) + im_true, im_test = _as_floats(im_true, im_test) + + norm_type = norm_type.lower() + if norm_type == 'euclidean': + denom = np.sqrt(np.mean((im_true*im_true), dtype=np.float64)) + elif norm_type == 'min-max': + denom = im_true.max() - im_true.min() + elif norm_type == 'mean': + denom = im_true.mean() + else: + raise ValueError("Unsupported norm_type") + return np.sqrt(mean_squared_error(im_true, im_test)) / denom + + +def psnr(im_true, im_test, dynamic_range=None): + """ Compute the peak signal to noise ratio (PSNR) for an image. + + Parameters + ---------- + im_true : ndarray + Ground-truth image. + im_test : ndarray + Test image. + dynamic_range : int + The dynamic range of the input image (distance between minimum and + maximum possible values). By default, this is estimated from the image + data-type. + + Returns + ------- + psnr : float + The PSNR metric. + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio + + """ + _assert_compatible(im_true, im_test) + if dynamic_range is None: + dmin, dmax = dtype_range[im_true.dtype.type] + true_min, true_max = np.min(im_true), np.max(im_true) + if true_max > dmax or true_min < dmin: + raise ValueError( + "im_true has intensity values outside the range expected for " + "its data type. Please manually specify the dynamic_range") + if true_min >= 0: + # most common case (255 for uint8, 1 for float) + dynamic_range = dmax + else: + dynamic_range = dmax - dmin + + im_true, im_test = _as_floats(im_true, im_test) + + err = mean_squared_error(im_true, im_test) + return 10 * np.log10((dynamic_range ** 2) / err) diff --git a/skimage/measure/tests/test_simple_metrics.py b/skimage/measure/tests/test_simple_metrics.py new file mode 100644 index 00000000..b59a11a8 --- /dev/null +++ b/skimage/measure/tests/test_simple_metrics.py @@ -0,0 +1,61 @@ +import numpy as np +from numpy.testing import (run_module_suite, assert_equal, assert_raises, + assert_almost_equal) + +from skimage.measure import psnr, normalized_root_mse, mean_squared_error +import skimage.data + +np.random.seed(5) +cam = skimage.data.camera() +sigma = 20.0 +cam_noisy = np.clip(cam + sigma * np.random.randn(*cam.shape), 0, 255) +cam_noisy = cam_noisy.astype(cam.dtype) + + +def test_PSNR_vs_IPOL(): + # Tests vs. imdiff result from the following IPOL article and code: + # http://www.ipol.im/pub/art/2011/g_lmii/ + p_IPOL = 22.4497 + p = psnr(cam, cam_noisy) + assert_almost_equal(p, p_IPOL, decimal=4) + + +def test_PSNR_float(): + p_uint8 = psnr(cam, cam_noisy) + p_float64 = psnr(cam/255., cam_noisy/255., dynamic_range=1) + assert_almost_equal(p_uint8, p_float64, decimal=5) + + +def test_PSNR_errors(): + assert_raises(ValueError, psnr, cam, cam.astype(np.float32)) + assert_raises(ValueError, psnr, cam, cam[:-1, :]) + + +def test_NRMSE(): + x = np.ones(4) + y = np.asarray([0., 2., 2., 2.]) + assert_equal(normalized_root_mse(y, x, 'mean'), 1/np.mean(y)) + assert_equal(normalized_root_mse(y, x, 'Euclidean'), 1/np.sqrt(3)) + assert_equal(normalized_root_mse(y, x, 'min-max'), 1/(y.max()-y.min())) + + +def test_NRMSE_no_int_overflow(): + camf = cam.astype(np.float32) + cam_noisyf = cam_noisy.astype(np.float32) + assert_almost_equal(mean_squared_error(cam, cam_noisy), + mean_squared_error(camf, cam_noisyf)) + assert_almost_equal(normalized_root_mse(cam, cam_noisy), + normalized_root_mse(camf, cam_noisyf)) + + +def test_NRMSE_errors(): + x = np.ones(4) + assert_raises(ValueError, normalized_root_mse, + x.astype(np.uint8), x.astype(np.float32)) + assert_raises(ValueError, normalized_root_mse, x[:-1], x) + # invalid normalization name + assert_raises(ValueError, normalized_root_mse, x, x, 'foo') + + +if __name__ == "__main__": + run_module_suite()