psnrhvsm.m

   1 function [p_hvs_m, p_hvs] = psnrhvsm(img1, img2, wstep)\r
   2 \r
   3 %========================================================================\r
   4 %\r
   5 % Calculation of PSNR-HVS-M and PSNR-HVS image quality measures\r
   6 %\r
   7 % PSNR-HVS-M is Peak Signal to Noise Ratio taking into account \r
   8 % Contrast Sensitivity Function (CSF) and between-coefficient   \r
   9 % contrast masking of DCT basis functions\r
  10 % PSNR-HVS is Peak Signal to Noise Ratio taking into account only CSF \r
  11 %\r
  12 % Copyright(c) 2006 Nikolay Ponomarenko \r
  13 % All Rights Reserved\r
  14 %\r
  15 % Homepage: http://ponomarenko.info, E-mail: nikolay{}ponomarenko.info\r
  16 %\r
  17 %----------------------------------------------------------------------\r
  18 %\r
  19 % Permission to use, copy, or modify this software and its documentation\r
  20 % for educational and research purposes only and without fee is hereby\r
  21 % granted, provided that this copyright notice and the original authors'\r
  22 % names appear on all copies and supporting documentation. This program\r
  23 % shall not be used, rewritten, or adapted as the basis of a commercial\r
  24 % software or hardware product without first obtaining permission of the\r
  25 % authors. The authors make no representations about the suitability of\r
  26 % this software for any purpose. It is provided "as is" without express\r
  27 % or implied warranty.\r
  28 %\r
  29 %----------------------------------------------------------------------\r
  30 %\r
  31 % This is an implementation of the algorithm for calculating the PSNR-HVS-M\r
  32 % or PSNR-HVS between two images. Please refer to the following papers:\r
  33 %\r
  34 % PSNR-HVS-M:\r
  35 % [1] Nikolay Ponomarenko, Flavia Silvestri, Karen Egiazarian, Marco Carli, \r
  36 %     Jaakko Astola, Vladimir Lukin, "On between-coefficient contrast masking \r
  37 %     of DCT basis functions", CD-ROM Proceedings of the Third International \r
  38 %     Workshop on Video Processing and Quality Metrics for Consumer Electronics \r
  39 %     VPQM-07, Scottsdale, Arizona, USA, 25-26 January, 2007, 4 p.\r
  40 %\r
  41 % PSNR-HVS:\r
  42 % [2] K. Egiazarian, J. Astola, N. Ponomarenko, V. Lukin, F. Battisti, \r
  43 %     M. Carli, New full-reference quality metrics based on HVS, CD-ROM \r
  44 %     Proceedings of the Second International Workshop on Video Processing \r
  45 %     and Quality Metrics, Scottsdale, USA, 2006, 4 p.\r
  46 %\r
  47 % Kindly report any suggestions or corrections to uagames{}mail.ru\r
  48 %\r
  49 %----------------------------------------------------------------------\r
  50 %\r
  51 % Input : (1) img1: the first image being compared\r
  52 %         (2) img2: the second image being compared\r
  53 %         (3) wstep: step of 8x8 window to calculate DCT \r
  54 %             coefficients. Default value is 8.\r
  55 %\r
  56 % Output: (1) p_hvs_m: the PSNR-HVS-M value between 2 images.\r
  57 %             If one of the images being compared is regarded as \r
  58 %             perfect quality, then PSNR-HVS-M can be considered as the\r
  59 %             quality measure of the other image.\r
  60 %             If compared images are visually undistingwished, \r
  61 %             then PSNR-HVS-M = 100000.\r
  62 %         (2) p_hvs: the PSNR-HVS value between 2 images.\r
  63 %\r
  64 % Default Usage:\r
  65 %   Given 2 test images img1 and img2, whose dynamic range is 0-255\r
  66 %\r
  67 %   [p_hvs_m, p_hvs] = psnrhvsm(img1, img2);\r
  68 %\r
  69 % See the results:\r
  70 %\r
  71 %   p_hvs_m  % Gives the PSNR-HVS-M value\r
  72 %   p_hvs    % Gives the PSNR-HVS value\r
  73 %\r
  74 %========================================================================\r
  75 \r
  76 if nargin < 2\r
  77   p_hvs_m = -Inf;\r
  78   p_hvs = -Inf;\r
  79   return;\r
  80 end\r
  81 \r
  82 if size(img1) ~= size(img2)\r
  83   p_hvs_m = -Inf;\r
  84   p_hvs = -Inf;\r
  85   return;\r
  86 end\r
  87 \r
  88 if nargin > 2 \r
  89   step = wstep;\r
  90 else\r
  91   step = 8; % Default value is 8;\r
  92 end\r
  93 \r
  94 img1=double(img1);\r
  95 img2=double(img2);\r
  96 \r
  97 LenXY=size(img1);LenY=LenXY(1);LenX=LenXY(2);\r
  98 \r
  99 CSFCof  = [1.608443, 2.339554, 2.573509, 1.608443, 1.072295, 0.643377, 0.504610, 0.421887;\r
 100            2.144591, 2.144591, 1.838221, 1.354478, 0.989811, 0.443708, 0.428918, 0.467911;\r
 101            1.838221, 1.979622, 1.608443, 1.072295, 0.643377, 0.451493, 0.372972, 0.459555;\r
 102            1.838221, 1.513829, 1.169777, 0.887417, 0.504610, 0.295806, 0.321689, 0.415082;\r
 103            1.429727, 1.169777, 0.695543, 0.459555, 0.378457, 0.236102, 0.249855, 0.334222;\r
 104            1.072295, 0.735288, 0.467911, 0.402111, 0.317717, 0.247453, 0.227744, 0.279729;\r
 105            0.525206, 0.402111, 0.329937, 0.295806, 0.249855, 0.212687, 0.214459, 0.254803;\r
 106            0.357432, 0.279729, 0.270896, 0.262603, 0.229778, 0.257351, 0.249855, 0.259950];\r
 107 % see an explanation in [2]\r
 108 \r
 109 MaskCof = [0.390625, 0.826446, 1.000000, 0.390625, 0.173611, 0.062500, 0.038447, 0.026874;\r
 110            0.694444, 0.694444, 0.510204, 0.277008, 0.147929, 0.029727, 0.027778, 0.033058;\r
 111            0.510204, 0.591716, 0.390625, 0.173611, 0.062500, 0.030779, 0.021004, 0.031888;\r
 112            0.510204, 0.346021, 0.206612, 0.118906, 0.038447, 0.013212, 0.015625, 0.026015;\r
 113            0.308642, 0.206612, 0.073046, 0.031888, 0.021626, 0.008417, 0.009426, 0.016866;\r
 114            0.173611, 0.081633, 0.033058, 0.024414, 0.015242, 0.009246, 0.007831, 0.011815;\r
 115            0.041649, 0.024414, 0.016437, 0.013212, 0.009426, 0.006830, 0.006944, 0.009803;\r
 116            0.019290, 0.011815, 0.011080, 0.010412, 0.007972, 0.010000, 0.009426, 0.010203];\r
 117 % see an explanation in [1]\r
 118 \r
 119 S1 = 0; S2 = 0; Num = 0;\r
 120 X=1;Y=1;\r
 121 while Y <= LenY-7\r
 122   while X <= LenX-7 \r
 123     A = img1(Y:Y+7,X:X+7);\r
 124     B = img2(Y:Y+7,X:X+7);\r
 125     A_dct = dct2(A); B_dct = dct2(B);\r
 126     MaskA = maskeff(A,A_dct);\r
 127     MaskB = maskeff(B,B_dct);\r
 128     if MaskB > MaskA\r
 129       MaskA = MaskB;\r
 130     end\r
 131     X = X + step;\r
 132     for k = 1:8\r
 133       for l = 1:8\r
 134         u = abs(A_dct(k,l)-B_dct(k,l));\r
 135         S2 = S2 + (u*CSFCof(k,l)).^2;    % PSNR-HVS\r
 136         if (k~=1) | (l~=1)               % See equation 3 in [1]\r
 137           if u < MaskA/MaskCof(k,l) \r
 138             u = 0;\r
 139           else\r
 140             u = u - MaskA/MaskCof(k,l);\r
 141           end\r
 142         end\r
 143         S1 = S1 + (u*CSFCof(k,l)).^2;    % PSNR-HVS-M\r
 144         Num = Num + 1;\r
 145       end\r
 146     end\r
 147   end\r
 148   X = 1; Y = Y + step;\r
 149 end\r
 150 \r
 151 if Num ~=0\r
 152   S1 = S1/Num;S2 = S2/Num;\r
 153   if S1 == 0 \r
 154     p_hvs_m = 100000; % img1 and img2 are visually undistingwished\r
 155   else\r
 156     p_hvs_m = 10*log10(255*255/S1);\r
 157   end\r
 158   if S2 == 0  \r
 159     p_hvs = 100000; % img1 and img2 are identical\r
 160   else\r
 161     p_hvs = 10*log10(255*255/S2);\r
 162   end\r
 163 end\r
 164 \r
 165 function m=maskeff(z,zdct)  \r
 166 % Calculation of Enorm value (see [1])\r
 167 m = 0;\r
 168 \r
 169 MaskCof = [0.390625, 0.826446, 1.000000, 0.390625, 0.173611, 0.062500, 0.038447, 0.026874;\r
 170            0.694444, 0.694444, 0.510204, 0.277008, 0.147929, 0.029727, 0.027778, 0.033058;\r
 171            0.510204, 0.591716, 0.390625, 0.173611, 0.062500, 0.030779, 0.021004, 0.031888;\r
 172            0.510204, 0.346021, 0.206612, 0.118906, 0.038447, 0.013212, 0.015625, 0.026015;\r
 173            0.308642, 0.206612, 0.073046, 0.031888, 0.021626, 0.008417, 0.009426, 0.016866;\r
 174            0.173611, 0.081633, 0.033058, 0.024414, 0.015242, 0.009246, 0.007831, 0.011815;\r
 175            0.041649, 0.024414, 0.016437, 0.013212, 0.009426, 0.006830, 0.006944, 0.009803;\r
 176            0.019290, 0.011815, 0.011080, 0.010412, 0.007972, 0.010000, 0.009426, 0.010203];\r
 177 % see an explanation in [1]\r
 178 \r
 179 for k = 1:8\r
 180   for l = 1:8\r
 181     if (k~=1) | (l~=1)\r
 182       m = m + (zdct(k,l).^2) * MaskCof(k,l);\r
 183     end\r
 184   end\r
 185 end\r
 186 pop=vari(z);\r
 187 if pop ~= 0\r
 188   pop=(vari(z(1:4,1:4))+vari(z(1:4,5:8))+vari(z(5:8,5:8))+vari(z(5:8,1:4)))/pop;\r
 189 end\r
 190 m = sqrt(m*pop)/32;   % sqrt(m*pop/16/64)\r
 191 \r
 192 function d=vari(AA)\r
 193   d=var(AA(:))*length(AA(:));\r