psnrhma.m

function [p_hvs_m, p_hvs] = psnrhma(img1, img2, wstep)

%========================================================================
%
% Calculation of PSNR-HMA and PSNR-HA image quality measures
%
% PSNR-HMA is Peak Signal to Noise Ratio taking into account
% Contrast Sensitivity Function (CSF), between-coefficient
% contrast masking of DCT basis functions as well as mean shift and
% contrast changing
% PSNR-HA is the same metric without taking into account constrast masking
%
% Copyright(c) 2011 Nikolay Ponomarenko
% All Rights Reserved
%
% Homepage: www.ponomarenko.info , E-mail: nikolay@ponomarenko.info
%
%----------------------------------------------------------------------
%
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereby
% granted, provided that this copyright notice and the original authors
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for calculating the PSNR-HA
% or PSNR-HMA between two images. Please refer to the following paper:
%
% N. Ponomarenko, O. Ieremeiev, V. Lukin, K. Egiazarian, M. Carli,
% Modified Image Visual Quality Metrics for Contrast Change and Mean Shift
% Accounting, Proceedings of CADSM, February 2011, Ukraine, pp. 305 - 311.
%
% Kindly report any suggestions or corrections to nikolay@ponomarenko.info
%
%----------------------------------------------------------------------
%
% Input : (1) img1: the first image being compared
%         (2) img2: the second image being compared
%         (3) wstep: step of 8x8 window to calculate DCT
%             coefficients. Default value is 8.
%
% Output: (1) p_hvs_m: the PSNR-HMA value between 2 images.
%             If one of the images being compared is regarded as
%             perfect quality, then PSNR-HMA can be considered as the
%             quality measure of the other image.
%             If compared images are visually undistingwished,
%             then PSNR-HMA = 100000.
%         (2) p_hvs: the PSNR-HA value between 2 images.
%
% Default Usage:
%   Given 2 test images img1 and img2, whose dynamic range is 0-255
%
%   [p_hvs_m, p_hvs] = psnrhvsm(img1, img2);
%
% See the results:
%
%   p_hvs_m  % Gives the PSNR-HMA value
%   p_hvs    % Gives the PSNR-HA value
%
%========================================================================

  if nargin < 2
    p_hvs_m = -Inf;
    p_hvs = -Inf;
    return;
  end

  if size(img1) ~= size(img2)
    p_hvs_m = -Inf;
    p_hvs = -Inf;
    return;
  end

  if nargin > 2
    step = wstep;
  else
    step = 8; % Default value is 8;
  end

  CSFCof_y  = [1.6193873005, 2.2901594831, 2.08509755623, 1.48366094411, 1.00227514334, 0.678296995242, 0.466224900598, 0.3265091542;
           2.2901594831, 1.94321815382, 2.04793073064, 1.68731108984, 1.2305666963, 0.868920337363, 0.61280991668, 0.436405793551;
           2.08509755623, 2.04793073064, 1.34329019223, 1.09205635862, 0.875748795257, 0.670882927016, 0.501731932449, 0.372504254596;
           1.48366094411, 1.68731108984, 1.09205635862, 0.772819797575, 0.605636379554, 0.48309405692, 0.380429446972, 0.295774038565;
           1.00227514334, 1.2305666963, 0.875748795257, 0.605636379554, 0.448996256676, 0.352889268808, 0.283006984131, 0.226951348204;
           0.678296995242, 0.868920337363, 0.670882927016, 0.48309405692, 0.352889268808, 0.27032073436, 0.215017739696, 0.17408067321;
           0.466224900598, 0.61280991668, 0.501731932449, 0.380429446972, 0.283006984131, 0.215017739696, 0.168869545842, 0.136153931001;
           0.3265091542, 0.436405793551, 0.372504254596, 0.295774038565, 0.226951348204, 0.17408067321, 0.136153931001, 0.109083846276];

  CSFCof_cb  = [1.91113096927, 2.46074210438, 1.18284184739, 1.14982565193, 1.05017074788, 0.898018824055, 0.74725392039, 0.615105596242;
           2.46074210438, 1.58529308355, 1.21363250036, 1.38190029285, 1.33100189972, 1.17428548929, 0.996404342439, 0.830890433625;
           1.18284184739, 1.21363250036, 0.978712413627, 1.02624506078, 1.03145147362, 0.960060382087, 0.849823426169, 0.731221236837;
           1.14982565193, 1.38190029285, 1.02624506078, 0.861317501629, 0.801821139099, 0.751437590932, 0.685398513368, 0.608694761374;
           1.05017074788, 1.33100189972, 1.03145147362, 0.801821139099, 0.676555426187, 0.605503172737, 0.55002013668, 0.495804539034;
           0.898018824055, 1.17428548929, 0.960060382087, 0.751437590932, 0.605503172737, 0.514674450957, 0.454353482512, 0.407050308965;
           0.74725392039, 0.996404342439, 0.849823426169, 0.685398513368, 0.55002013668, 0.454353482512, 0.389234902883, 0.342353999733;
           0.615105596242, 0.830890433625, 0.731221236837, 0.608694761374, 0.495804539034, 0.407050308965, 0.342353999733, 0.295530605237];

  CSFCof_cr  = [2.03871978502, 2.62502345193, 1.26180942886, 1.11019789803, 1.01397751469, 0.867069376285, 0.721500455585, 0.593906509971;
           2.62502345193, 1.69112867013, 1.17180569821, 1.3342742857, 1.28513006198, 1.13381474809, 0.962064122248, 0.802254508198;
           1.26180942886, 1.17180569821, 0.944981930573, 0.990876405848, 0.995903384143, 0.926972725286, 0.820534991409, 0.706020324706;
           1.11019789803, 1.3342742857, 0.990876405848, 0.831632933426, 0.77418706195, 0.725539939514, 0.661776842059, 0.587716619023;
           1.01397751469, 1.28513006198, 0.995903384143, 0.77418706195, 0.653238524286, 0.584635025748, 0.531064164893, 0.478717061273;
           0.867069376285, 1.13381474809, 0.926972725286, 0.725539939514, 0.584635025748, 0.496936637883, 0.438694579826, 0.393021669543;
           0.721500455585, 0.962064122248, 0.820534991409, 0.661776842059, 0.531064164893, 0.438694579826, 0.375820256136, 0.330555063063;
           0.593906509971, 0.802254508198, 0.706020324706, 0.587716619023, 0.478717061273, 0.393021669543, 0.330555063063, 0.285345396658];

  [y,x,z]=size(img1);
  if z>1
    a=rgb2ycbcr(img1);
    a1=double(a(:,:,1));
    a2=double(a(:,:,2));
    a3=double(a(:,:,3));
    b=rgb2ycbcr(img2);
    b1=double(b(:,:,1));
    b2=double(b(:,:,2));
    b3=double(b(:,:,3));

    [p11,p12]=msehma(a1, b1, step, CSFCof_y, 0, 0.37);
    [p21,p22]=msehma(a2, b2, step, CSFCof_cb, 0.25, 2);
    [p31,p32]=msehma(a3, b3, step, CSFCof_cr, 0.25, 2);

    S1=(p11+(p21+p31)*0.5)/2;
    S2=(p12+(p22+p32)*0.5)/2;
  else
    [S1, S2]=msehma(img1, img2, step, CSFCof_y, 0, 0.37);
  end

  if S1 == 0
    p_hvs_m = 100000; % img1 and img2 are visually undistingwished
  else
    p_hvs_m = 10*log10(255*255/S1);
  end
  if S2 == 0
    p_hvs = 100000; % img1 and img2 are identical
  else
    p_hvs = 10*log10(255*255/S2);
  end

function [S1, S2] = msehma(img1, img2, step, CSFCof, KofContr1, KofContr2)

  LenXY=size(img1);LenY=LenXY(1);LenX=LenXY(2);

  MaskCof = (CSFCof.*0.3885746225901003).^2;
  % see an explanation in [1]

  delt=(sum(double(img1(:)))-sum(double(img2(:))))/length(img1(:));
  img2m=double(img2)+delt;

  mean1=mean(img1(:));
  mean2=mean(img2m(:));
  tmp=(img1-mean1).*(img2m-mean2);
  sq=var(img2m(:),1);
  l=sum(tmp(:))/length(tmp(:))/sq;
  if l<1
    KofContr=KofContr1;
  else
    KofContr=KofContr2;
  end;
  img3m=mean2+(img2m-mean2)*l;
  S3=(sum(sum((img2m-img1).^2))-sum(sum((img3m-img1).^2)))/length(img1(:));

  S1 = 0; S2 = 0; Num = 0;
  SS1=0;SS2=0;
  X=1;Y=1;
  while Y <= LenY-7
    while X <= LenX-7
      A = img1(Y:Y+7,X:X+7);
      B = img2m(Y:Y+7,X:X+7);
      B2 = img3m(Y:Y+7,X:X+7);
      A_dct = dct2(A); B_dct = dct2(B); B_dct2 = dct2(B2);
      MaskA = maskeff(A,A_dct,MaskCof);
      MaskB = maskeff(B,B_dct,MaskCof);
      if MaskB > MaskA
        MaskA = MaskB;
      end
      X = X + step;
      for k = 1:8
        for l = 1:8
          u = abs(A_dct(k,l)-B_dct(k,l));
          u2= abs(A_dct(k,l)-B_dct2(k,l));
          S2 = S2 + (u*CSFCof(k,l)).^2;    % PSNR-HVS
          SS2 = SS2 + (u2*CSFCof(k,l)).^2;    % PSNR-HVS

          if (k~=1) || (l~=1)               % See equation 3 in [1]
            if u < MaskA/MaskCof(k,l)
              u = 0;
            else
              u = u - MaskA/MaskCof(k,l);
            end;
            if u2 < MaskA/MaskCof(k,l)
              u2 = 0;
            else
              u2 = u2 - MaskA/MaskCof(k,l);
            end;
          end
          S1 = S1 + (u*CSFCof(k,l)).^2;    % PSNR-HVS-M
          SS1 = SS1 + (u2*CSFCof(k,l)).^2;    % PSNR-HVS-M
          Num = Num + 1;
        end
      end
    end
    X = 1; Y = Y + step;
  end

  if Num ~=0
    S1 = S1/Num;S2 = S2/Num;
    SS1 = SS1/Num;SS2=SS2/Num;
    delt=delt^2;

    if (S1>SS1) S1=SS1+(S1-SS1)*KofContr;
    end
    S1=S1 + 0.045*delt;
    if (S2>SS2) S2=SS2+(S2-SS2)*KofContr;
    end
    S2=S2 + 0.045*delt;
  end

function m=maskeff(z,zdct,MaskCof)
% Calculation of Enorm value (see [1])
  m = 0;
  for k = 1:8
    for l = 1:8
      if (k~=1) || (l~=1)
        m = m + (zdct(k,l).^2) * MaskCof(k,l);
      end
    end
  end
  pop=vari(z);
  if pop ~= 0
    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;
  end
  m = sqrt(m*pop)/32;   % sqrt(m*pop/16/64)

function d=vari(AA)
  d=var(AA(:))*length(AA(:));