function ratio = wsnr(orig, dith, nfreq) %WSNR Weighted signal to noise ratio. % RATIO = WSNR(IM1, IM2, MAXFREQ) computes the weighted signal to noise % ratio of IM2 with respect to IM1 and returns the result in dB. % A weighting appropriate to the human visual system, as proposed % by Mannos and Sakrison (1) and modified by Mitsa and Varkur (2), % is used. MAXFREQ specifies the spatial frequency in cycles per % degree that corresponds to the Nyquist frequency in the x-direction. % If MAXFREQ is omitted, it defaults to 60 cyc/deg. % % Refs: (1) J. Mannos and D. Sakrison, "The effects of a visual % fidelity criterion on the encoding of images", IEEE Trans. Inf. % Theory, IT-20(4), pp. 525-535, July 1974. % % (2) T. Mitsa and K. Varkur, "Evaluation of contrast sensitivity % functions for the formulation of quality measures incorporated % in halftoning algorithms", IEEE Int. Conference on Acoustics, Speech % and Signal Processing '93-V, pp. 301-304. % % See also SNR, PSNR, TSNR, SNRVSF. % From the Image Quality Assessment Toolbox, N. Damera-Venkata and % Prof. Brian L. Evans, released April 28, 2001 if nargin==2 nfreq=60; end [x,y]=size(orig); [xplane,yplane]=meshgrid(-x/2+0.5:x/2-0.5); % generate mesh plane=(xplane+i*yplane)/x*2*nfreq; radfreq=abs(plane); % radial frequency % According to (2), we modify the radial frequency according to angle. % w is a symmetry parameter that gives approx. 3 dB down along the % diagonals. w=0.7; s=(1-w)/2*cos(4*angle(plane))+(1+w)/2; radfreq=radfreq./s; % Now generate the CSF csf=2.6*(0.0192+0.114*radfreq).*exp(-(0.114*radfreq).^1.1); f=find(radfreq<7.8909); csf(f)=0.9809+zeros(size(f)); % Find weighted SNR in frequency domain. Note that, because we are not % weighting the signal, we compute signal power in the spatial domain. % This requires us to multiply by the image size in pixels to get the % signal power in the frequency domain for division. err=orig-dith; %err_wt=fft2(err).*csf; % weighted error spectrum err_wt=fftshift(fft2(err)).*csf; im=fft2(orig); mse=sum(sum(err_wt.*conj(err_wt))); % weighted error power mss=sum(sum(im.*conj(im))); % signal power ratio=10*log10(mss/mse); % compute SNR