Merge pull request #1897 from grlee77/simple_metrics

FEAT: Simple image comparison metrics (PSNR, NRMSE)

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']

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)

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()