diff --git a/pytorch/__pycache__/examplesLDR.cpython-38.pyc b/pytorch/__pycache__/examplesLDR.cpython-38.pyc new file mode 100644 index 0000000..0e74770 Binary files /dev/null and b/pytorch/__pycache__/examplesLDR.cpython-38.pyc differ diff --git a/pytorch/__pycache__/pu21_display_model.cpython-38.pyc b/pytorch/__pycache__/pu21_display_model.cpython-38.pyc new file mode 100644 index 0000000..a223d22 Binary files /dev/null and b/pytorch/__pycache__/pu21_display_model.cpython-38.pyc differ diff --git a/pytorch/__pycache__/pu21_encoder.cpython-38.pyc b/pytorch/__pycache__/pu21_encoder.cpython-38.pyc new file mode 100644 index 0000000..4a2f360 Binary files /dev/null and b/pytorch/__pycache__/pu21_encoder.cpython-38.pyc differ diff --git a/pytorch/__pycache__/pu21_metric.cpython-38.pyc b/pytorch/__pycache__/pu21_metric.cpython-38.pyc new file mode 100644 index 0000000..832e94f Binary files /dev/null and b/pytorch/__pycache__/pu21_metric.cpython-38.pyc differ diff --git a/pytorch/examples.py b/pytorch/examples.py new file mode 100644 index 0000000..42ff2c4 --- /dev/null +++ b/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) ) diff --git a/pytorch/examplesLDR.py b/pytorch/examplesLDR.py new file mode 100644 index 0000000..251b4a9 --- /dev/null +++ b/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) ) diff --git a/pytorch/psnr.py b/pytorch/psnr.py new file mode 100644 index 0000000..13f0376 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/pytorch/pu21_display_model.py b/pytorch/pu21_display_model.py new file mode 100644 index 0000000..5a2d611 --- /dev/null +++ b/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 )) \ No newline at end of file diff --git a/pytorch/pu21_encoder.py b/pytorch/pu21_encoder.py new file mode 100644 index 0000000..49b1104 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pytorch/pu21_metric.py b/pytorch/pu21_metric.py new file mode 100644 index 0000000..0666216 --- /dev/null +++ b/pytorch/pu21_metric.py @@ -0,0 +1,90 @@ +import torch +from pu21_encoder import pu21_encoder +from ignite.metrics import SSIM,PSNR +from skimage.metrics import structural_similarity,peak_signal_noise_ratio + +def pu21_metric( I_test, I_reference, metric, display_model=None ): + """ + A convenience function for calling traditional (SDR) metrics on + PU-encoded pixel values. This is useful for adapting traditional metrics + to HDR images. + + Q = pu21_metric( I_test, I_reference, metric ) % for HDR images + Q = pu21_metric( I_test, I_reference, metric, display_model ) % for SDR + + When no display model is passed, I_test and I_reference must be provided + as ABSOLUTE linear colour or luminance values. If unsure what those values + are, please refer to the paper [1]. + + When display model is passed, I_test and I_reference contain images in + display-encoded sRGB colour space or luma (standard images). If images + are stored as floating point values, the values must be in the range 0-1. + If they are stored as integers, they must be in the range + from 0 to maxint(class(I)). + + display_model is an object of the class pu21_display_model_gog. Check + exaples/ex_sdr_images.m on how to create an object of this class. + + metric can be either 'PSNR', 'SSIM'. You can also pass a handle to a + function Q = fun(I_test, I_referece). The function should expect both + images to be stored as floating point numbers and the typical range from + 0 to 255. The range of values can be higher for bright HDR images or + bright SDR displays. + + [1] 1. Mantiuk RK. + Practicalities of predicting quality of high dynamic range images and video. + In: 2016 IEEE International Conference on Image Processing (ICIP), p. 904–8. + https://www.cl.cam.ac.uk/~rkm38/pdfs/mantiuk2016prac_hdr_metrics.pdf + http://dx.doi.org/10.1109/ICIP.2016.7532488 + + + If images are stored as intigers, convert to a single precision floating + point between 0 and 1. + """ + I_test = torch.tensor(I_test) + if not I_test.is_floating_point(): + print( Warning("Hardcoded 255, this is super cursed to do in python")) + I_test = I_test.to(torch.float64)/torch.iinfo(I_test.dtype).max + if not I_reference.is_floating_point(): + print( Warning("Hardcoded 255, this is super cursed to do in python")) + print(I_reference.dtype) + I_reference = I_reference.to(torch.float64)/torch.iinfo(I_reference.dtype).max + # I_test/=2**16-1 + # I_reference/=2**16-1 + I_reference = I_reference.to(torch.float64) + I_test = I_test.to(torch.float64) + if display_model !=None: + L_test = display_model.forward( I_test ) + L_reference = display_model.forward( I_reference ) + else: + L_test = I_test + L_reference = I_reference + + pu21 = pu21_encoder() + P_test = pu21.encode( L_test ) + P_reference = pu21.encode( L_reference ) + P_test = P_test.permute(2,1,0) + P_reference = P_reference.permute(2,1,0) + print("P_test max", torch.max( P_test)) + print("P_test min", torch.min( P_test)) + if isinstance(metric, str): + metricFunc = None + if metric.lower() == 'psnr': + #metricFunc = PSNR(1.0) + return peak_signal_noise_ratio(P_test.numpy(),P_reference.numpy(),data_range=1) + if metric.lower() == 'ssim': + # Note that we are passing floating point values, which are in the + # range 0-256 for the luminance range 0.1 to 100 cd/m^2 + #metricFunc = SSIM(255,kernel_size=(101,101)) + P_test = P_test.permute(1,2,0) + P_reference = P_reference.permute(1,2,0) + return structural_similarity(P_test.numpy(),P_reference.numpy(),multichannel=True,gaussian_weights=True, use_sample_covariance=False,data_range=255.0) + if metricFunc==None: + raise Exception( 'Unknown metric {}'.format(metric) ) + else: + metricFunc.update((P_test.unsqueeze(0),P_reference.unsqueeze(0))) + Q = metricFunc.compute() + else: + Q = metric(P_test,P_reference) + return Q +