Merge pull request #119 from stefanv/ssim

ENH: Add structural similarity index and gradient.

doc/examples/plot_ssim.py

@@ -0,0 +1,68 @@
+'''
+===========================
+Structural similarity index
+===========================
+
+When comparing images, the mean squared error (MSE)--while simple to
+implement--is not highly indicative of perceived similarity.  Structural
+similarity aims to address this shortcoming by taking texture into account
+[1]_, [2]_.
+
+The example shows two modifications of the input image, each with the same MSE,
+but with very different mean structural similarity indices.
+
+.. [1] Zhou Wang; Bovik, A.C.; ,"Mean squared error: Love it or leave it? A new
+       look at Signal Fidelity Measures," Signal Processing Magazine, IEEE,
+       vol. 26, no. 1, pp. 98-117, Jan. 2009.
+
+.. [2] Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, "Image quality
+       assessment: From error visibility to structural similarity," IEEE
+       Transactions on Image Processing, vol. 13, no. 4, pp. 600-612,
+       Apr. 2004.
+
+'''
+
+from skimage import data, color, io, exposure, img_as_float
+from skimage.measure import structural_similarity as ssim
+
+import numpy as np
+
+img = img_as_float(data.camera())
+rows, cols = img.shape
+
+noise = np.ones_like(img) * 0.2 * (img.max() - img.min())
+noise[np.random.random(size=noise.shape) > 0.5] *= -1
+
+def mse(x, y):
+    return np.linalg.norm(x - y)
+
+img_noise = img + noise
+img_const = img + abs(noise)
+
+import matplotlib.pyplot as plt
+f, (ax0, ax1, ax2) = plt.subplots(1, 3)
+
+mse_none = mse(img, img)
+ssim_none = ssim(img, img, dynamic_range=img.max() - img.min())
+
+mse_noise = mse(img, img_noise)
+ssim_noise = ssim(img, img_noise, dynamic_range=img_const.max() - img_const.min())
+
+mse_const = mse(img, img_const)
+ssim_const = ssim(img, img_const, dynamic_range=img_noise.max() - img_noise.min())
+
+label = 'MSE: %2.f, SSIM: %.2f'
+
+ax0.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=1)
+ax0.set_xlabel(label % (mse_none, ssim_none))
+ax0.set_title('Original image')
+
+ax1.imshow(img_noise, cmap=plt.cm.gray, vmin=0, vmax=1)
+ax1.set_xlabel(label % (mse_noise, ssim_noise))
+ax1.set_title('Image with noise')
+
+ax2.imshow(img_const, cmap=plt.cm.gray, vmin=0, vmax=1)
+ax2.set_xlabel(label % (mse_const, ssim_const))
+ax2.set_title('Image plus constant')
+
+plt.show()

skimage/measure/__init__.py

@@ -1,2 +1,3 @@
 from .find_contours import find_contours
 from ._regionprops import regionprops
+from ._structural_similarity import structural_similarity

skimage/measure/_structural_similarity.py

@@ -0,0 +1,105 @@
+from __future__ import division
+
+__all__ = ['structural_similarity']
+
+import numpy as np
+
+from ..util.dtype import dtype_range
+from ..util.shape import view_as_windows
+
+def structural_similarity(X, Y, win_size=7,
+                          gradient=False, dynamic_range=None):
+    """Compute the mean structural similarity index between two images.
+
+    Parameters
+    ----------
+    X, Y : (N,N) ndarray
+        Images.
+    win_size : int
+        The side-length of the sliding window used in comparison.  Must
+        be an odd value.
+    gradient : bool
+        If True, also return the gradient.
+    dynamic_range : int
+        Dynamic range of the input image (distance between minimum and
+        maximum possible values).  By default, this is estimated from
+        the image data-type.
+
+    Returns
+    -------
+    s : float
+        Strucutural similarity.
+    grad : (N * N,) ndarray
+        Gradient of the structural similarity index between X and Y.
+        This is only returned if `gradient` is set to True.
+
+    References
+    ----------
+    .. [1] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P.
+       (2004). Image quality assessment: From error visibility to
+       structural similarity. IEEE Transactions on Image Processing,
+       13, 600-612.
+
+    """
+    if not X.dtype == Y.dtype:
+        raise ValueError('Input images must have the same dtype.')
+
+    if not X.shape == Y.shape:
+        raise ValueError('Input images must have the same dimensions.')
+
+    if not (win_size % 2 == 1):
+        raise ValueError('Window size must be odd.')
+
+    if dynamic_range is None:
+        dmin, dmax = dtype_range[X.dtype.type]
+        dynamic_range = dmax - dmin
+
+    XW = view_as_windows(X, (win_size, win_size))
+    YW = view_as_windows(Y, (win_size, win_size))
+
+    NS = len(XW)
+    NP = win_size * win_size
+
+    ux = np.mean(np.mean(XW, axis=2), axis=2)
+    uy = np.mean(np.mean(YW, axis=2), axis=2)
+
+    # Compute variances var(X), var(Y) and var(X, Y)
+    cov_norm = 1 / (win_size**2 - 1)
+    XWM = XW - ux[..., None, None]
+    YWM = YW - uy[..., None, None]
+    vx = cov_norm * np.sum(np.sum(XWM**2, axis=2), axis=2)
+    vy = cov_norm * np.sum(np.sum(YWM**2, axis=2), axis=2)
+    vxy = cov_norm * np.sum(np.sum(XWM * YWM, axis=2), axis=2)
+
+    R = dynamic_range
+    K1 = 0.01
+    K2 = 0.03
+    C1 = (K1 * R)**2
+    C2 = (K2 * R)**2
+
+    A1, A2, B1, B2 = (v[..., None, None] for v in
+                      (2 * ux * uy + C1,
+                       2 * vxy + C2,
+                       ux**2 + uy**2 + C1,
+                       vx + vy + C2))
+
+    S = np.mean((A1 * A2) / (B1 * B2))
+
+    if gradient:
+        local_grad = 2 / (NP * B1**2 * B2**2) * \
+            (
+            A1 * B1 * (B2 * XW - A2 * YW) - \
+            B1 * B2 * (A2 - A1) * ux[..., None, None] + \
+            A1 * A2 * (B1 - B2) * uy[..., None, None]
+            )
+
+        grad = np.zeros_like(X, dtype=float)
+        OW = view_as_windows(grad, (win_size, win_size))
+
+        OW += local_grad
+        grad /= NS
+
+        return S, grad
+
+    else:
+        return S

skimage/measure/tests/test_structural_similarity.py

@@ -0,0 +1,58 @@
+import numpy as np
+from numpy.testing import assert_equal
+
+from skimage.measure import structural_similarity as ssim
+import scipy.optimize as opt
+
+def test_ssim_patch_range():
+    N = 51
+    X = (np.random.random((N, N)) * 255).astype(np.uint8)
+    Y = (np.random.random((N, N)) * 255).astype(np.uint8)
+
+    assert(ssim(X, Y, win_size=N) < 0.1)
+    assert_equal(ssim(X, X, win_size=N), 1)
+
+def test_ssim_image():
+    N = 100
+    X = (np.random.random((N, N)) * 255).astype(np.uint8)
+    Y = (np.random.random((N, N)) * 255).astype(np.uint8)
+
+    S0 = ssim(X, X, win_size=3)
+    assert_equal(S0, 1)
+
+    S1 = ssim(X, Y, win_size=3)
+    assert(S1 < 0.3)
+
+## Come up with a better way of testing the gradient
+##
+## def test_ssim_grad():
+##     N = 30
+##     X = np.random.random((N, N)) * 255
+##     Y = np.random.random((N, N)) * 255
+
+##     def func(Y):
+##         return ssim(X, Y, dynamic_range=255)
+
+##     def grad(Y):
+##         return ssim(X, Y, dynamic_range=255, gradient=True)[1]
+
+##     assert(np.all(opt.check_grad(func, grad, Y) < 0.05))
+
+def test_ssim_dtype():
+    N = 30
+    X = np.random.random((N, N))
+    Y = np.random.random((N, N))
+
+    S1 = ssim(X, Y)
+
+    X = (X * 255).astype(np.uint8)
+    Y = (X * 255).astype(np.uint8)
+
+    S2 = ssim(X, Y)
+
+    assert S1 < 0.1
+    assert S2 < 0.1
+
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()

skimage/util/__init__.py

@@ -1 +1,2 @@
 from .dtype import *
+from .shape import *