Merge pull request #13 from up42/metrics_review

Metrics review

Author: markusuwe

README.md

@@ -2,16 +2,14 @@
 
 Implementation of eight evaluation metrics to access the similarity between two images. The eight metrics are as follows:
 
-<i><a href="https://en.wikipedia.org/wiki/Root-mean-square_deviation">Root mean square error (RMSE)</a></i>,
-<i><a href="https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio">Peak signal-to-noise ratio (PSNR)</a></i>,
-<i><a href="https://en.wikipedia.org/wiki/Structural_similarity">Structural Similarity Index (SSIM)</a></i>,
-<i><a href="https://www.tandfonline.com/doi/full/10.1080/22797254.2019.1628617">Information theoretic-based Statistic Similarity Measure (ISSM)</a></i>,
-<i><a href="https://www4.comp.polyu.edu.hk/~cslzhang/IQA/TIP_IQA_FSIM.pdf">Feature-based similarity index (FSIM)</a></i>,
-<i><a href="https://www.sciencedirect.com/science/article/abs/pii/S0924271618302636">Signal to reconstruction error ratio (SRE)</a></i>,
-<i><a href="https://ntrs.nasa.gov/citations/19940012238">Spectral angle mapper (SAM)</a></i>, and
-<i><a href="https://www.researchgate.net/publication/3342733_A_Universal_Image_Quality_Index">Universal image quality index (UIQ)</a></i>
-
-
+ * <i><a href="https://en.wikipedia.org/wiki/Root-mean-square_deviation">Root mean square error (RMSE)</a></i>,
+ * <i><a href="https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio">Peak signal-to-noise ratio (PSNR)</a></i>,
+ * <i><a href="https://en.wikipedia.org/wiki/Structural_similarity">Structural Similarity Index (SSIM)</a></i>,
+ * <i><a href="https://www4.comp.polyu.edu.hk/~cslzhang/IQA/TIP_IQA_FSIM.pdf">Feature-based similarity index (FSIM)</a></i>,
+ * <i><a href="https://www.tandfonline.com/doi/full/10.1080/22797254.2019.1628617">Information theoretic-based Statistic Similarity Measure (ISSM)</a></i>,
+ * <i><a href="https://www.sciencedirect.com/science/article/abs/pii/S0924271618302636">Signal to reconstruction error ratio (SRE)</a></i>,
+ * <i><a href="https://ntrs.nasa.gov/citations/19940012238">Spectral angle mapper (SAM)</a></i>, and
+ * <i><a href="https://ece.uwaterloo.ca/~z70wang/publications/quality_2c.pdf">Universal image quality index (UIQ)</a></i>
 
 ## Instructions
 

image_similarity_measures/evaluate.py

@@ -78,6 +78,7 @@ def evaluation(org_img_path, pred_img_path, mode, metric, write_to_file):
         metric_dict[metric] = {f"{metric.upper()}": out_value}
         write_final_dict(metric, metric_dict)
 
+
 def main():
     parser = argparse.ArgumentParser(description="Evaluates an Image Super Resolution Model")
     parser.add_argument("--org_img_path", type=str, help="Path to original input image")

image_similarity_measures/quality_metrics.py

@@ -17,37 +17,40 @@ def _assert_image_shapes_equal(org_img: np.ndarray, pred_img: np.ndarray, metric
     assert org_img.shape == pred_img.shape, msg
 
 
-def rmse(org_img: np.ndarray, pred_img: np.ndarray, data_range=4096):
+def rmse(org_img: np.ndarray, pred_img: np.ndarray, max_p=4095) -> float:
     """
     Root Mean Squared Error
-    """
 
+    Calculated individually for all bands, then averaged
+    """
     _assert_image_shapes_equal(org_img, pred_img, "RMSE")
-    rmse_final = []
+
+    rmse_bands = []
     for i in range(org_img.shape[2]):
-        m = np.mean(np.square((org_img[:, :, i] - pred_img[:, :, i]) / data_range))
+        m = np.mean(np.square((org_img[:, :, i] - pred_img[:, :, i]) / max_p))
         s = np.sqrt(m)
-        rmse_final.append(s)
-    return np.mean(rmse_final)
+        rmse_bands.append(s)
 
+    return np.mean(rmse_bands)
 
-def psnr(org_img: np.ndarray, pred_img: np.ndarray, data_range=4096):
+
+def psnr(org_img: np.ndarray, pred_img: np.ndarray, max_p=4095) -> float:
     """
-    Peek Signal to Noise Ratio, a measure similar to mean squared error.
+    Peek Signal to Noise Ratio, implemented as mean squared error converted to dB.
 
     It can be calculated as
     PSNR = 20 * log10(MAXp) - 10 * log10(MSE)
 
-    When using 12-bit imagery MaxP is 4096, for 8-bit imagery 256
+    When using 12-bit imagery MaxP is 4095, for 8-bit imagery 255. For floating point imagery using values between
+    0 and 1 (e.g. unscaled reflectance) the first logarithmic term can be dropped as it becomes 0
     """
     _assert_image_shapes_equal(org_img, pred_img, "PSNR")
 
-    r = []
+    mse_bands = []
     for i in range(org_img.shape[2]):
-        val = 20 * np.log10(data_range) - 10. * np.log10(np.mean(np.square(org_img[:, :, i] - pred_img[:, :, i])))
-        r.append(val)
+        mse_bands.append(np.mean(np.square(org_img[:, :, i] - pred_img[:, :, i])))
 
-    return np.mean(r)
+    return 20 * np.log10(max_p) - 10. * np.log10(np.mean(mse_bands))
 
 
 def _similarity_measure(x, y, constant):
@@ -70,7 +73,7 @@ def _gradient_magnitude(img: np.ndarray, img_depth):
     return np.sqrt(scharrx ** 2 + scharry ** 2)
 
 
-def fsim(org_img: np.ndarray, pred_img: np.ndarray):
+def fsim(org_img: np.ndarray, pred_img: np.ndarray, T1=0.85, T2=160) -> float:
     """
     Feature-based similarity index, based on phase congruency (PC) and image gradient magnitude (GM)
 
@@ -80,11 +83,22 @@ def fsim(org_img: np.ndarray, pred_img: np.ndarray):
 
     There are also alternatives to implement GM, the FSIM authors suggest to use the Scharr
     operation which is implemented in OpenCV.
+
+    Note that FSIM is defined in the original papers for grayscale as well as for RGB images. Our use cases
+    are mostly multi-band images e.g. RGB + NIR. To accommodate for this fact, we compute FSIM for each individual
+    band and then take the average.
+
+    Note also that T1 and T2 are constants depending on the dynamic range of PC/GM values. In theory this parameters
+    would benefit from fine-tuning based on the used data, we use the values found in the original paper as defaults.
+
+    Args:
+        org_img -- numpy array containing the original image
+        pred_img -- predicted image
+        T1 -- constant based on the dynamic range of PC values
+        T2 -- constant based on the dynamic range of GM values
     """
     _assert_image_shapes_equal(org_img, pred_img, "FSIM")
 
-    T1 = 0.85  # a constant based on the dynamic range PC
-    T2 = 160  # a constant based on the dynamic range GM
     alpha = beta = 1  # parameters used to adjust the relative importance of PC and GM features
     fsim_list = []
     for i in range(org_img.shape[2]):
@@ -145,9 +159,12 @@ def _edge_c(x, y):
     return numerator / denominator
 
 
-def issm(org_img: np.ndarray, pred_img: np.ndarray):
+def issm(org_img: np.ndarray, pred_img: np.ndarray) -> float:
     """
     Information theoretic-based Statistic Similarity Measure
+
+    Note that the term e which is added to both the numerator as well as the denominator is not properly
+    introduced in the paper. We assume the authers refer to the Euler number.
     """
     _assert_image_shapes_equal(org_img, pred_img, "ISSM")
 
@@ -167,13 +184,13 @@ def issm(org_img: np.ndarray, pred_img: np.ndarray):
     return np.nan_to_num(numerator / denominator)
 
 
-def ssim(org_img: np.ndarray, pred_img: np.ndarray, data_range=4096):
+def ssim(org_img: np.ndarray, pred_img: np.ndarray, max_p=4095) -> float:
     """
     Structural SIMularity index
     """
     _assert_image_shapes_equal(org_img, pred_img, "SSIM")
 
-    return structural_similarity(org_img, pred_img, data_range=data_range, multichannel=True)
+    return structural_similarity(org_img, pred_img, data_range=max_p, multichannel=True)
 
 
 def sliding_window(image, stepSize, windowSize):
@@ -184,49 +201,58 @@ def sliding_window(image, stepSize, windowSize):
             yield (x, y, image[y:y + windowSize[1], x:x + windowSize[0]])
 
 
-def uiq(org_img: np.ndarray, pred_img: np.ndarray):
+def uiq(org_img: np.ndarray, pred_img: np.ndarray, step_size=1, window_size=8):
     """
     Universal Image Quality index
     """
     # TODO: Apply optimization, right now it is very slow
     _assert_image_shapes_equal(org_img, pred_img, "UIQ")
     q_all = []
-    for (x, y, window_org), (x, y, window_pred) in zip(sliding_window(org_img, stepSize=1, windowSize=(8, 8)),
-                                                       sliding_window(pred_img, stepSize=1, windowSize=(8, 8))):
+    for (x, y, window_org), (x, y, window_pred) in zip(sliding_window(org_img, stepSize=step_size,
+                                                                      windowSize=(window_size, window_size)),
+                                                       sliding_window(pred_img, stepSize=step_size,
+                                                                      windowSize=(window_size, window_size))):
         # if the window does not meet our desired window size, ignore it
         if window_org.shape[0] != 8 or window_org.shape[1] != 8:
             continue
-        org_img_mean = np.mean(org_img)
-        pred_img_mean = np.mean(pred_img)
-        org_img_variance = np.var(org_img)
-        pred_img_variance = np.var(pred_img)
-        org_pred_img_variance = np.mean((window_org - org_img_mean) * (window_pred - pred_img_mean))
 
-        numerator = 4 * org_pred_img_variance * org_img_mean * pred_img_mean
-        denominator = (org_img_variance + pred_img_variance) * (org_img_mean**2 + pred_img_mean**2)
+        for i in range(org_img.shape[2]):
+            org_band = window_org[:, :, i]
+            pred_band = window_pred[:, :, i]
+            org_band_mean = np.mean(org_band)
+            pred_band_mean = np.mean(pred_band)
+            org_band_variance = np.var(org_band)
+            pred_band_variance = np.var(pred_band)
+            org_pred_band_variance = np.mean((org_band - org_band_mean) * (pred_band - pred_band_mean))
 
-        if denominator != 0.0:
-            q = numerator / denominator
-            q_all.append(q)
+            numerator = 4 * org_pred_band_variance * org_band_mean * pred_band_mean
+            denominator = (org_band_variance + pred_band_variance) * (org_band_mean**2 + pred_band_mean**2)
+
+            if denominator != 0.0:
+                q = numerator / denominator
+                q_all.append(q)
 
     return np.mean(q_all)
 
 
-def sam(org_img: np.ndarray, pred_img: np.ndarray):
+def sam(org_img: np.ndarray, pred_img: np.ndarray, convert_to_degree=True):
     """
-    calculates spectral angle mapper
+    Spectral Angle Mapper which defines the spectral similarity between two spectra
     """
+
     _assert_image_shapes_equal(org_img, pred_img, "SAM")
-    org_img = org_img.reshape((org_img.shape[0] * org_img.shape[1], org_img.shape[2]))
-    pred_img = pred_img.reshape((pred_img.shape[0] * pred_img.shape[1], pred_img.shape[2]))
 
-    N = org_img.shape[1]
-    sam_angles = np.zeros(N)
-    for i in range(org_img.shape[1]):
-        val = np.clip(np.dot(org_img[:, i], pred_img[:, i]) / (np.linalg.norm(org_img[:, i]) * np.linalg.norm(pred_img[:, i])), -1, 1)
-        sam_angles[i] = np.arccos(val)
+    # Spectral angles are first computed for each pair of pixels
+    numerator = np.sum(np.multiply(pred_img, org_img), axis=2)
+    denominator = np.linalg.norm(org_img, axis=2) * np.linalg.norm(pred_img, axis=2)
+    val = np.clip(numerator / denominator, -1, 1)
+    sam_angles = np.arccos(val)
+    if convert_to_degree:
+        sam_angles = sam_angles * 180.0 / np.pi
 
-    return np.mean(sam_angles * 180.0 / np.pi)
+    # The original paper states that SAM values are expressed as radians, while e.g. Lanares
+    # et al. (2018) use degrees. We therefore made this configurable, with degree the default
+    return np.mean(sam_angles)
 
 
 def sre(org_img: np.ndarray, pred_img: np.ndarray):
@@ -238,8 +264,8 @@ def sre(org_img: np.ndarray, pred_img: np.ndarray):
     sre_final = []
     for i in range(org_img.shape[2]):
         numerator = np.square(np.mean(org_img[:, :, i]))
-        denominator = ((np.linalg.norm(org_img[:, :, i] - pred_img[:, :, i]))) /\
+        denominator = (np.linalg.norm(org_img[:, :, i] - pred_img[:, :, i])) /\
                       (org_img.shape[0] * org_img.shape[1])
-        sre_final.append(10 * np.log10(numerator/denominator))
+        sre_final.append(numerator/denominator)
 
-    return np.mean(sre_final)
+    return 10 * np.log10(np.mean(sre_final))

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="image-similarity-measures",
-    version="0.2.2",
+    version="0.3.0",
     author="UP42",
     author_email="[email protected]",
     description="Evaluation metrics to assess the similarity between two images.",