functionf = WPSNR(A,B,varargin) % ThisfunctioncomputesWPSNR(weightedpeaksignal - to - noiseratio)between % twoimages.Theanswer is in decibels(dB). % % Usingcontrastsensitivityfunction(CSF)toweightspatialfrequency % oferrorimage. % % Using:WPSNR(A,B) % % WrittenbyRuizhenLiu,http: // www.assuredigit.com if A == Berror( ' Imagesareidentical:PSNRhasinfinitevaluelilizong ' )endmax2_A = max(max(A));max2_B = max(max(B));min2_A = min(min(A));min2_B = min(min(B)); if max2_A > 1 | max2_B > 1 | min2_A < 0 | min2_B < 0 error( ' inputmatricesmusthavevaluesintheinterval[0,1] ' )ende = A - B; if nargin < 3 fc = csf; % filtercoefficientsofCSF else fc = varargin ... {1} ;endew = filter2(fc,e); % filteringerrorwithCSFdecibels = 20 * log10( 1 / (sqrt(mean(mean(ew. ^ 2 ))))); % disp(sprintf( ' WPSNR=+%5.2fdB ' ,decibels))f = decibels; %============= functionfc = csf() %============= % ProgramtocomputeCSF % ComputecontrastsensitivityfunctionofHVS % % Output:fc --- filtercoefficientsofCSF % % Reference: % MakotoMiyahara % " ObjectivePictureQualityScale(PQS)forImageCoding " % IEEETrans.onComm.,Vol 46 ,No. 9 , 1998 . % % WrittenbyRuizhenLiu,http: // www.assuredigit.com % computefrequencyresponsematrixFmat = csfmat; % Plotfrequencyresponse % mesh(Fmat);pause % compute 2 - Dfiltercoefficient using FSAMP2fc = fsamp2(Fmat); % mesh(fc) %=================== functionFmat = csfmat() %=================== % ComputeCSFfrequencyresponsematrix % % Calledbyfunctioncsf.m % % Frequencyrange:therangoffrequencyseemstobe: % w = pi = ( 2 * pi * f) / 60 % f = 60 * w / ( 2 * pi),about 21.2 % min_f = - 20 ;max_f = 20 ;step_f = 1 ;u = min_f:step_f:max_f;v = min_f:step_f:max_f;n = length(u);Z = zeros(n); for i = 1 :n for j = 1 :nZ(i,j) = csffun(u(i),v(j)); % callingfunctioncsffunendendFmat = Z; %======================== functionSa = csffun(u,v) %======================== % ContrastSensitivityFunction in spatialfrequency % Thisfilecomputethespatialfrequencyweightingoferrors % % Reference: % MakotoMiyahara % " ObjectivePictureQualityScale(PQS)forImageCoding " % IEEETrans.onComm.,Vol 46 ,No. 9 , 1998 . % % Input:u --- horizontalspatialfrequencies % v --- verticalspatialfrequencies % % Output:frequencyresponse % % WrittenbyRuizhenLiu,http: // www.assuredigit.com % ComputeSa -- spatialfrequencyresponse % symsSwsigmafuvsigma = 2 ;f = sqrt(u. * u + v. * v);w = 2 * pi * f / 60 ;Sw = 1.5 * exp( - sigma ^ 2 * w ^ 2 / 2 ) - exp( - 2 * sigma ^ 2 * w ^ 2 / 2 ); % Modification in Highfrequencysita = atan(v. / (u + eps));bita = 8 ;f0 = 11.13 ;w0 = 2 * pi * f0 / 60 ;Ow = ( 1 + exp(bita * (w - w0)) * (cos( 2 * sita)) ^ 4 ) / ( 1 + exp(bita * (w - w0))); % ComputefinalresponseSa = Sw * Ow;