From 3529e4d8188e890b0988a3f7674092e08c12e4fb Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Tue, 8 Nov 2011 17:27:22 -0800 Subject: [PATCH] ENH: Add Structural SIMilarity (SSIM) image comparison. --- skimage/measure/__init__.py | 1 + skimage/measure/_ssim.py | 110 +++++++++++++++++++++++++++++ skimage/measure/tests/test_ssim.py | 36 ++++++++++ 3 files changed, 147 insertions(+) create mode 100644 skimage/measure/_ssim.py create mode 100644 skimage/measure/tests/test_ssim.py diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 422db569..cb9dc679 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -1,2 +1,3 @@ from .find_contours import find_contours from ._regionprops import regionprops +from ._ssim import ssim diff --git a/skimage/measure/_ssim.py b/skimage/measure/_ssim.py new file mode 100644 index 00000000..a93571e8 --- /dev/null +++ b/skimage/measure/_ssim.py @@ -0,0 +1,110 @@ +from __future__ import division + +__all__ = ['ssim'] + +import numpy as np +from numpy.lib import stride_tricks + +def _as_windows(X, win_size=7): + """Re-stride an array to simulate a sliding window. + + Parameters + ---------- + X : 2D-ndarray + Input image. + + Returns + ------- + window : (N, win_size, win_size) ndarray + Sliding windows. + + """ + if not X.ndim == 2: + raise ValueError('Input images must be 2-dimensional.') + + X = np.ascontiguousarray(X) + r, c = X.shape + + strides = X.strides + row_jump, el_jump = strides + half_width = (win_size // 2) + + new_strides = (row_jump, el_jump, row_jump, el_jump) + new_rows = r - 2 * half_width + new_cols = c - 2 * half_width + new_shape = (new_rows, new_cols, win_size, win_size) + + windows = stride_tricks.as_strided(X, shape=new_shape, strides=new_strides) + windows = windows.reshape((-1, win_size, win_size)) + + return windows + + +def ssim(X, Y, win_size=7, dynamic_range=255): + """Compute the 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. + dynamic_range : int + Dynamic range of the input image (distance between minimum and + maximum possible values). This should eventually be + auto-computed, but just specifying it manually for now. + + Returns + ------- + s : float + Strucutural similarity. + + 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('Inout images must have the same dimensions.') + + import time + + tic = time.time() + + XW = _as_windows(X, win_size=win_size) + YW = _as_windows(Y, win_size=win_size) + + tic = time.time() + + # Flatten windows + XW = XW.reshape(XW.shape[0], -1) + YW = YW.reshape(YW.shape[0], -1) + + ux = np.mean(XW, axis=1) + uy = np.mean(YW, axis=1) + + tic = time.time() + + # Compute variances var(X), var(Y) and var(X, Y) + cov_norm = 1 / (win_size**2 - 1) + XWM = XW - ux[:, None] + YWM = YW - uy[:, None] + vx = cov_norm * np.sum(XWM**2, axis=1) + vy = cov_norm * np.sum(YWM**2, axis=1) + vxy = cov_norm * np.sum(XWM * YWM, axis=1) + + R = dynamic_range + K1 = 0.01 + K2 = 0.03 + C1 = (K1 * R)**2 + C2 = (K2 * R)**2 + + return np.mean(((2 * ux * uy + C1) * (2 * vxy + C2)) / \ + ((ux**2 + uy**2 + C1) * (vx + vy + C2))) diff --git a/skimage/measure/tests/test_ssim.py b/skimage/measure/tests/test_ssim.py new file mode 100644 index 00000000..8e6684b0 --- /dev/null +++ b/skimage/measure/tests/test_ssim.py @@ -0,0 +1,36 @@ +import numpy as np +from numpy.testing import assert_equal + +from skimage.measure._ssim import ssim, _as_windows + +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_as_windows(): + X = np.arange(100).reshape((10, 10)) + W = _as_windows(X, win_size=7) + assert_equal(len(W), 16) + + W = _as_windows(X, win_size=3) + assert_equal(W[0], [[0, 1, 2], + [10, 11, 12], + [20, 21, 22]]) + +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) + +if __name__ == "__main__": + np.testing.run_module_suite()