% This script compares the WInHD inverse halftoning estimate with
% Multi-scale gradient estimate of Kite et al.
%
% Reference: 
% "WInHD: Wavelet-based Inverse Halftoning via Deconvolution",
% R. Neelamani, R. Nowak, R. Baraniuk, Submitted to the IEEE
% Transactions on Image Processing, September 2002. Available for
% download from http://www.dsp.rice.edu/publications/
%


% Inverse halftoning of Floyd halftones

RunMe;

x = double(imread('lena.tif')); % Original Lena image
y = double(imread('lena_floyd_ht.tif')); % Floyd-Stienberg halftone

Grad_estimate = double(imread('lena_floyd_Grad_iht.tif')); 

% Gradestimate obtained from 
% http://www.ece.utexas.edu/~bevans/projects/inverseHalftoning.html
%
% \bibitem{gradient_halftoning}
%  T.~D. Kite, N.~Damera-Venkata, B.~L. Evans, and A.~C. Bovik, ``A high 
% quality, fast inverse halftoning algorithm for error diffused 
% halftones, '' {\em Proc. IEEE Int. Conf. Image Processing -- ICIP~'98}
% , vol.~2, pp.~59--63, Oct. 4--7, 1998.

% WInHD_estimate obtained by inverse halftoning y with the IHT_Floyd
% function, but with complex wavelets (not included with software)
WInHD_estimate_cwt = double(imread('lena_floyd_cwt_WInHD_iht.tif')); 

% WInHD_estimate obtained by inverse halftoning y with the IHT_Floyd function.
WInHD_estimate = double(imread('lena_floyd_WInHD_iht.tif')); 


psnr_winhd_cwt = psnr(WInHD_estimate_cwt,x);
psnr_winhd = psnr(WInHD_estimate,x);
psnr_grad = psnr(Grad_estimate,x);

figure(1)
clf
subplot(2,2,1)
imagesc(x)
imset
im_range = get(gca,'CLim');
title('Original')

subplot(2,2,2)
imagesc(y)
imset
title('Halftoned  with Floyd')
imset

subplot(2,2,3)
imagesc(Grad_estimate,im_range)
imset
title(strcat('Gradient estimate,','PSNR = ',num2str(psnr_grad),'dB'))

subplot(2,2,4)
imagesc(WInHD_estimate_cwt,im_range)
imset
title(strcat('WInHD estimate (complex wavelet based),','PSNR = ',num2str(psnr_winhd_cwt),'dB'))

figure(2)
clf
subplot(2,2,1)
imagesc(x)
imset
im_range = get(gca,'CLim');
title('Original')

subplot(2,2,2)
imagesc(y)
imset
title('Halftoned  with Floyd')
imset

subplot(2,2,3)
imagesc(Grad_estimate,im_range)
imset
title(strcat('Gradient estimate,','PSNR = ',num2str(psnr_grad),'dB'))

subplot(2,2,4)
imagesc(WInHD_estimate,im_range)
imset
title(strcat('WInHD estimate (redundant wavelet based),','PSNR = ',num2str(psnr_winhd),'dB'))

figure(1)

disp(['To compare visual qualities, please download the Halftoning ToolBox  from http://www.ece.utexas.edu/~bevans/projects/halftoning/index.html. The routines wsnr.m and img_qi.m are used to report the WSNR and Universal Image Quality values in the Journal paper.'])

% WInHD_WSNR = wsnr(WInHD_estimate_cwt,x,60)
% Grad_WSNR = wsnr(Grad_estimate,x,60)
% WInHD_UIQI = img_qi(WInHD_estimate_cwt,x)
% Grad_UIQI = img_qi(Grad_estimate,x)

