%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the results of the WaRD algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
subplot(2,2,2)
imagesc(X)
colormap('gray')
im_range = get(gca,'CLim');
imagesc(Y, im_range)
axis('square')
rmaxis
title(strcat('Noisy image, SNR = ', num2str(BSNR),'dB'))

subplot(2,2,1)
imagesc(X, im_range)
axis('square')
rmaxis
title(strcat('Original image'))

subplot(2,2,3)
imagesc(Wiener_Estimate, im_range)
axis('square')
rmaxis
ISNR_Wiener = 10*log10(norm(X(:)-Y(:))^2 /norm(Wiener_Estimate(:) - X(:))^2 );
SNR_Wiener = 10*log10(norm(X(:))^2 /norm(Wiener_Estimate(:) - X(:))^2 );
title(strcat('Wiener, ISNR =',num2str(ISNR_Wiener),'dB, SNR =',num2str(SNR_Wiener),'dB'))

subplot(2,2,4)
imagesc(x_ward, im_range)
axis('square')
rmaxis
ISNR_WDWF = 10*log10(norm(X(:) - Y(:))^2 /norm(x_ward(:) - X(:))^2 );
SNR_WDWF = 10*log10(norm(X(:))^2 /norm(x_ward(:) - X(:))^2 );
title(strcat('WaRD, ISNR = ', num2str(ISNR_WDWF), 'dB, SNR = ', num2str(SNR_WDWF),'dB'))

