sample vs. population covariance difference addressed to more accurately match reference implementations

skimage/measure/_structural_similarity.py

@@ -77,7 +77,8 @@ def _discard_edges(X, pad):
 def structural_similarity(X, Y, win_size=None, gradient=False,
                           dynamic_range=None, multichannel=None,
                           gaussian_weights=False, full=False,
-                          image_content_weighting=False):
+                          image_content_weighting=False,
+                          use_sample_covariance=True):
     """Compute the mean structural similarity index between two images.
 
     Parameters
@@ -106,6 +107,9 @@ def structural_similarity(X, Y, win_size=None, gradient=False,
     image_content_weighting : bool
         If True, weight the ssim mean is spatially weighted by image content as
         proposed in Wang and Shang 2006 [3].
+    use_sample_covariance : bool
+        if True, normalize covariances by N-1 rather than, N where N is the
+        number of pixels within the sliding window.
 
     Returns
     -------
@@ -120,7 +124,8 @@ def structural_similarity(X, Y, win_size=None, gradient=False,
     Notes
     -----
     To exactly match the implementation of Wang et. al. [1], set
-    `gaussian_weights` to True and `win_size` to 11.
+    `gaussian_weights` to True, `win_size` to 11, and `use_sample_covariance`
+    to False.
 
     References
     ----------
@@ -216,7 +221,12 @@ def structural_similarity(X, Y, win_size=None, gradient=False,
     Y = Y.astype(np.float64)
 
     NP = win_size ** ndim
-    cov_norm = NP / (NP - 1)  # filter will normalize by NP, but we want NP - 1
+
+    # filter has already normalized by NP
+    if use_sample_covariance:
+        cov_norm = NP / (NP - 1)  # sample covariance
+    else:
+        cov_norm = 1.0  # population covariance to match Wang et. al. 2004
 
     # compute (weighted) means
     ux = filter_func(X, **filter_args)

skimage/measure/tests/test_structural_similarity.py

@@ -105,8 +105,25 @@ def test_gaussian_mssim_vs_IPOL():
     # Tests vs. imdiff result from the following IPOL article and code:
     # http://www.ipol.im/pub/art/2011/g_lmii/
     mssim_IPOL = 0.327309966087341
-    mssim = ssim(cam, cam_noisy, gaussian_weights=True)
-    assert_almost_equal(mssim, mssim_IPOL, decimal=2)
+    mssim = ssim(cam, cam_noisy, gaussian_weights=True,
+                 use_sample_covariance=False)
+    assert_almost_equal(mssim, mssim_IPOL, decimal=5)
+
+
+def test_gaussian_mssim_vs_author_ref():
+    """
+    test vs. result from original author's Matlab implementation available at
+    https://ece.uwaterloo.ca/~z70wang/research/ssim/
+
+    Matlab test code:
+       img1 = imread('camera.png')
+       img2 = imread('camera_noisy.png')
+       mssim = ssim_index(img1, img2)
+    """
+    mssim_matlab = 0.218987555561590
+    mssim = ssim(cam, cam_noisy, gaussian_weights=True,
+                 use_sample_covariance=False)
+    assert_almost_equal(mssim, mssim_matlab, decimal=7)
 
 
 def test_gaussian_mssim_and_gradient_vs_Matlab():
@@ -118,9 +135,10 @@ def test_gaussian_mssim_and_gradient_vs_Matlab():
     grad_matlab = ref['grad_matlab']
     mssim_matlab = float(ref['mssim_matlab'])
 
-    mssim, grad = ssim(cam, cam_noisy, gaussian_weights=True, gradient=True)
+    mssim, grad = ssim(cam, cam_noisy, gaussian_weights=True, gradient=True,
+                       use_sample_covariance=False)
 
-    assert_almost_equal(mssim, mssim_matlab, decimal=2)
+    assert_almost_equal(mssim, mssim_matlab, decimal=7)
 
     # check almost equal aside from object borders
     assert_array_almost_equal(grad_matlab[5:-5], grad[5:-5])