Pytorch

pytorch/examples.py

@@ -0,0 +1,55 @@
+# This example shows how to run PU21 metrics on HDR images
+import HDRutils
+import os 
+import pu21_metric
+import torch
+import scipy.ndimage as ndimage
+import numpy as np
+import cv2
+# in python
+
+I_ref = HDRutils.imread( os.path.join("..","matlab","examples",'nancy_church.hdr' ))
+I_ref = torch.tensor(I_ref)
+print(I_ref.mean())
+L_peak = 4000 # Peak luminance of an HDR display
+
+# HDR images are often given in relative photometric units. They MUST be
+# mapped to absolute amount of light emitted from the display. For that, 
+# we map the peak value in the image to the peak value of the display:
+
+I_ref = I_ref/torch.max(I_ref) * L_peak
+
+# Add Gaussian noise of 20% contrast. Make sure all values are greater than
+# 0.05.
+
+I_test_noise = torch.maximum(I_ref + I_ref*torch.randn(I_ref.shape)*0.2, torch.tensor(0.05) )
+
+
+
+PSNR_noise = pu21_metric.pu21_metric( I_test_noise, I_ref, 'PSNR' )
+SSIM_noise = pu21_metric.pu21_metric( I_test_noise, I_ref, 'SSIM' )
+
+print('Image with noise: PSNR = {} dB, SSIM = {}'.format( PSNR_noise, SSIM_noise) )
+print(I_ref.shape)
+
+def matlab_style_gauss2D(shape=(3,3),sigma=0.5):
+    """
+    2D gaussian mask - should give the same result as MATLAB's
+    fspecial('gaussian',[shape],[sigma])
+    """
+    m,n = [(ss-1.)/2. for ss in shape]
+    y,x = np.ogrid[-m:m+1,-n:n+1]
+    h = np.exp( -(x*x + y*y) / (2.*sigma*sigma) )
+    h[ h < np.finfo(h.dtype).eps*h.max() ] = 0
+    sumh = h.sum()
+    if sumh != 0:
+        h /= sumh
+    return h
+I_test_blur = cv2.GaussianBlur(I_ref.numpy(), ksize=(0, 0), sigmaX=3, borderType=cv2.BORDER_REPLICATE)
+ #= torch.conv2d(I_ref.permute(2,1,0),torch.tensor(matlab_style_gauss2D()))
+#I_test_blur = ndimage.gaussian_filter(I_ref, sigma=3, order=0)
+#I_test_blur = imgaussfilt( I_ref, 3 )
+
+PSNR_noise = pu21_metric.pu21_metric( I_test_blur, I_ref, 'PSNR' )
+SSIM_noise = pu21_metric.pu21_metric( I_test_blur, I_ref, 'SSIM' )
+print('Image with blur: PSNR = {} dB, SSIM = {}'.format( PSNR_noise, SSIM_noise) )

pytorch/examplesLDR.py

@@ -0,0 +1,47 @@
+# This example shows how to run PU21 metrics on HDR images
+import HDRutils
+import os 
+import pu21_metric
+import pu21_display_model
+import torch
+import numpy as np
+import scipy.ndimage as ndimage
+import skimage
+import cv2
+
+I_ref = HDRutils.imread( os.path.join("..","matlab","examples",'wavy_facade.png' ))
+# maxVal = np.iinfo(np.int16).max
+I_ref = torch.tensor(I_ref.astype(np.int64))
+# print(I_ref.mean())
+# I_ref /=maxVal
+# print(I_ref.mean())
+
+# HDR images are often given in relative photometric units. They MUST be
+# mapped to absolute amount of light emitted from the display. For that, 
+# we map the peak value in the image to the peak value of the display:
+
+
+# Add Gaussian noise of 20% contrast. Make sure all values are greater than
+# 0.05.
+
+I_test_noise = torch.maximum( I_ref + I_ref*torch.randn(I_ref.shape)*0.001,torch.tensor(0))
+
+I_test_blur = cv2.GaussianBlur(I_ref.numpy().astype(np.float64), ksize=(0, 0), sigmaX=2, borderType=cv2.BORDER_REPLICATE)
+I_test_blur = I_test_blur / np.iinfo(np.uint16).max
+I_test_noise = I_test_noise / np.iinfo(np.uint16).max
+print("noised ",I_test_blur.mean())
+#I_test_blur = skimage.gaussian_filter(I_ref, sigma=2,mode = 'nearest',truncate=2.0)
+Y_peak = 100
+contrast = 1000
+E_ambient = 10
+pu_dm = pu21_display_model.pu21_display_model_gog( Y_peak, contrast, 2.2, E_ambient )
+
+PSNR_noise = pu21_metric.pu21_metric( I_test_noise, I_ref, 'PSNR',pu_dm )
+SSIM_noise = pu21_metric.pu21_metric( I_test_noise, I_ref, 'SSIM',pu_dm )
+
+
+PSNR_blur = pu21_metric.pu21_metric( I_test_blur, I_ref, 'PSNR',pu_dm )
+SSIM_blur = pu21_metric.pu21_metric( I_test_blur, I_ref, 'SSIM',pu_dm )
+
+print('Image with noise: PSNR = {} dB, SSIM = {}'.format( PSNR_noise, SSIM_noise) )
+print('Image with blur: PSNR = {} dB, SSIM = {}'.format( PSNR_blur, SSIM_blur) )

pytorch/psnr.py

@@ -0,0 +1,73 @@
+"""
+This code is from https://github.com/bonlime/pytorch-tools/blob/master/pytorch_tools/metrics/psnr.py
+"""
+
+import os
+import sys
+import math
+import torch
+import numpy as np
+import cv2
+
+
+class PSNR:
+    """Peak Signal to Noise Ratio
+    img1 and img2 have range [0, 255]"""
+
+    def __init__(self):
+        self.name = "PSNR"
+
+    @staticmethod
+    def __call__(img1, img2):
+        mse = torch.mean((img1 - img2) ** 2)
+        return 20 * torch.log10(255.0 / torch.sqrt(mse))
+
+
+class SSIM:
+    """Structure Similarity
+    img1, img2: [0, 255]"""
+
+    def __init__(self):
+        self.name = "SSIM"
+
+    @staticmethod
+    def __call__(self,img1, img2):
+        if not img1.shape == img2.shape:
+            raise ValueError("Input images must have the same dimensions.")
+        if img1.ndim == 2:  # Grey or Y-channel image
+            return self._ssim(img1, img2)
+        elif img1.ndim == 3:
+            if img1.shape[2] == 3:
+                ssims = []
+                for i in range(3):
+                    ssims.append(SSIM.ssim(img1, img2))
+                return np.array(ssims).mean()
+            elif img1.shape[2] == 1:
+                return self._ssim(np.squeeze(img1), np.squeeze(img2))
+        else:
+            raise ValueError("Wrong input image dimensions.")
+
+    @staticmethod
+    def _ssim(img1, img2):
+        C1 = (0.01 * 255) ** 2
+        C2 = (0.03 * 255) ** 2
+        img1 = img1.numpy()
+        img2 = img2.numpy()
+        img1 = img1.astype(np.float64)
+        img2 = img2.astype(np.float64)
+        kernel = cv2.getGaussianKernel(11, 1.5)
+        window = np.outer(kernel, kernel.transpose())
+
+        mu1 = cv2.filter2D(img1, -1, window)[5:-5, 5:-5]  # valid
+        mu2 = cv2.filter2D(img2, -1, window)[5:-5, 5:-5]
+        mu1_sq = mu1 ** 2
+        mu2_sq = mu2 ** 2
+        mu1_mu2 = mu1 * mu2
+        sigma1_sq = cv2.filter2D(img1 ** 2, -1, window)[5:-5, 5:-5] - mu1_sq
+        sigma2_sq = cv2.filter2D(img2 ** 2, -1, window)[5:-5, 5:-5] - mu2_sq
+        sigma12 = cv2.filter2D(img1 * img2, -1, window)[5:-5, 5:-5] - mu1_mu2
+
+        ssim_map = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2)) / (
+            (mu1_sq + mu2_sq + C1) * (sigma1_sq + sigma2_sq + C2)
+        )
+        return ssim_map.mean()

pytorch/pu21_display_model.py

@@ -0,0 +1,57 @@
+#%%
+import torch
+pi = torch.acos(torch.zeros(1)).item() * 2
+#%%
+class pu21_display_model_gog():
+    def __init__(self, Y_peak, contrast = 1000, gamma= 2.2, E_ambient= 0, k_refl = 0.005):
+        #        Gain-gamma-offset display model to simulate SDR displays
+
+        # dm = fvvdp_display_photo_gog( Y_peak )
+        # dm = fvvdp_display_photo_gog( Y_peak, contrast )
+        # dm = fvvdp_display_photo_gog( Y_peak, contrast, gamma )
+        # dm = fvvdp_display_photo_gog( Y_peak, contrast, gamma, E_ambient )
+        # dm = fvvdp_display_photo_gog( Y_peak, contrast, gamma, E_ambient, k_refl )
+
+        # Parameters (default value shown in []):
+        # Y_peak - display peak luminance in cd/m^2 (nit), e.g. 200 for a typical
+        #          office monitor
+        # contrast - [1000] the contrast of the display. The value 1000 means
+        #          1000:1
+        # gamma - [2.2] gamma of the display.
+        # E_ambient - [0] ambient light illuminance in lux, e.g. 600 for bright
+        #         office
+        # k_refl - [0.005] reflectivity of the display screen
+
+        # For more details on the GOG display model, see:
+        # https://www.cl.cam.ac.uk/~rkm38/pdfs/mantiuk2016perceptual_display.pdf
+
+        # Copyright (c) 2010-2021, Rafal Mantiuk 
+        self.Y_peak = Y_peak
+        self.contrast = contrast
+        self.gamma = gamma
+        self.E_ambient = E_ambient
+        self.k_refl = k_refl
+
+    def get_black_level(self):
+        Y_refl = self.E_ambient/pi*self.k_refl # Reflected ambient light            
+        Y_black = Y_refl + self.Y_peak/self.contrast
+        return Y_black
+
+    def forward(self,V):
+        # Transforms gamma-correctec pixel values V, which must be in the range
+        # 0-1, into absolute linear colorimetric values emitted from
+        # the display.
+        if torch.any(V>1) or torch.any(V<0) :
+            print( Warning('Pixel values must be in the range 0-1'))
+        Y_black = self.get_black_level()
+        L = (self.Y_peak-Y_black)*(V**self.gamma) + Y_black
+        return L
+
+    def print(self):
+        Y_black = self.get_black_level()
+        print( 'Photometric display model:' )
+        print( '  Peak luminance: {} cd/m^2'.format(self.Y_peak) )
+        print( '  Contrast - theoretical: {}:1'.format(round(self.contrast) ))
+        print( '  Contrast - effective: {}:1'.format(round(self.Y_peak/Y_black) ))
+        print( '  Ambient light: {} lux'.format(self.E_ambient ))
+        print( '  Display reflectivity: {}'.format(self.k_refl*100 ))

pytorch/pu21_encoder.py

@@ -0,0 +1,86 @@
+import torch
+
+class pu21_encoder():
+    """
+    Transform absolute linear luminance values to/from the perceptually
+    uniform (PU) space. This class is intended for adapting image quality 
+    metrics to operator on the HDR content.
+
+    Refer to the examples folder for examples how to use PU21 encoding
+    with both HDR and SDR/LDR images. 
+
+    The derivation of the PU21 encoding is explained in the paper: 
+
+    R. Mantiuk and M. Azimi
+    PU21: A novel perceptually uniform encoding for adapting existing
+    quality metrics for HDR.
+    Picture Coding Symposium 2021
+
+    The new PU21 encoding improves on the older PU (or PU08) encoding,
+    explained in: 
+
+    Aydin TO, Mantiuk R, Seidel H-P. 
+    Extending quality metrics to full luminance range images. 
+    In: Human Vision and Electronic Imaging. Spie 2008. no. 68060B. 
+    DOI: 10.1117/12.765095
+    """
+
+    def __init__(self,type = 'banding_glare'):
+        """ 
+        Create pu21_encoder or a given type (string), if the
+        parameter is supplied.
+
+        pu21 = pu21_encoder()
+        pu21 = pu21_encoder( type )
+
+        It is recommended that you use default type ('banding_glare')
+        by skipping 'type' parameter. 
+        """
+        self.L_min = 0.005
+        self.L_max = 10000
+        self.type = type
+        self.par = None
+        if type == 'banding':
+            self.par = [0,1.070275272, 0.4088273932, 0.153224308, 0.2520326168, 1.063512885, 1.14115047, 521.4527484]
+        if type == 'banding_glare':
+            self.par = [0,0.353487901, 0.3734658629, 8.277049286e-05, 0.9062562627, 0.09150303166, 0.9099517204, 596.3148142]
+        if type == 'peaks':
+            self.par = [0,1.043882782, 0.6459495343, 0.3194584211, 0.374025247, 1.114783422, 1.095360363, 384.9217577]
+        if type == 'peaks_glare':
+            self.par = [0,816.885024, 1479.463946, 0.001253215609, 0.9329636822, 0.06746643971, 1.573435413, 419.6006374]
+        if self.par ==None:
+            raise Exception("Unknow type: {}".format(type))
+
+    def encode(self,Y):
+        """
+        Convert from linear (optical) values Y to encoded (electronic) values V
+
+        V = encode(obj, Y)
+
+        V is in the range from 0 to 1.
+        Y is in the range from 0.005 to 10000. The values MUST be
+            scaled in the absolute units (nits, cd/m^2).
+        """
+        epsilon = 1e-5
+        if torch.any(Y<(self.L_min-epsilon)) or torch.any(Y>(self.L_max+epsilon)):
+            print( 'Values passed to encode are outside the valid range' )
+
+
+        Y = torch.clamp( Y, self.L_min, self.L_max); # Clamp the values
+        p = self.par
+        V = p[7] * (((p[1] + p[2]*(Y**p[4]) )/ (1+p[3]*(Y**p[4])))**p[5]-p[6])
+        return V
+
+    def decode(self, V):
+        """
+        Convert from encoded (electronic) values V into linear (optical) values Y
+
+        Y = decode(obj, V)
+
+        V is in the range from 0 to 1.
+        Y is in the range from 0.005 to 10000
+        """       
+        p = self.par
+        V_p = torch.maximum( V/p[7]+p[6], torch.tensor(0) )**(1/p[5])
+        Y = (torch.maximum( V_p-p[1], torch.tensor(0) ) /(p[2]-p[3]*V_p))**(1/p[4])
+        return Y