% Created by Ramesh Neelamani, Rice University on 11/9/99
%
% This file contains various examples to demonstrates the usage of the
% files in this package. 
%
% the experiments for peppers image can be run by substituting lena with
% peppers in the script below. the peppers image need to be downloaded 
% separately. 
%
% 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/
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


RunMe; % This simply adds the mtools and images directory to the current
       % matlab path. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 1 : Inverse Halftone Floyd halftoned image (Lena image)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Example 1: Inverse Halftone Floyd halftoned image (Lena image)')
x = double(imread('lena.tif')); 
y = double(imread('lena_floyd_ht.tif'));  % precomputed halftone
disp('Computing inverse halftone')
WInHD_estimate = IHT_Floyd(y,x);
disp('Pausing: Hit return to continue to next example')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 2 : Inverse Halftone Jarvis halftoned image (Lena image)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Example 2: Inverse Halftone Jarvis halftoned image (Lena image)')
x = double(imread('lena.tif')); 
y = double(imread('lena_jarvis_ht.tif'));  % precomputed halftone
disp('Computing inverse halftone')
WInHD_estimate = IHT_Jarvis(y,x);
disp('Pausing: Hit return to continue to next example')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 3: Computing Floyd halftone and inverse halftone (Lena image)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


disp('Example 3: Computing Floyd halftone and inverse halftone (Lena image)')
x = double(imread('lena.tif')); 
disp('Computing Floyd halftone')
y = floyd(x); 
disp('Computing inverse halftone')
WInHD_estimate = IHT_Floyd(y,x);
disp('Pausing: Hit return to continue to next example')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 4: Computing Jarvis halftone and inverse halftone (Lena image)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Example 2: Computing Jarvis halftone and inverse halftone (Lena image)')
x = double(imread('lena.tif')); 
disp('Computing Jarvis halftone')
y = jarvis(x); 
disp('Computing inverse halftone')
WInHD_estimate = IHT_Jarvis(y,x);
disp('Pausing: Hit return to continue to next example')
pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 5 : Inverse Halftone Floyd halftoned image (Peppers image)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Example 1: Inverse Halftone Floyd halftoned image (Peppers image)')
x = double(imread('peppers.tif')); 
y = double(imread('peppers_floyd_ht.tif'));  % precomputed halftone
disp('Computing inverse halftone')
WInHD_estimate = IHT_Floyd(y,x);
disp('Pausing: Hit return to continue to next example')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 5: Comparing the algorithm performance: See Compare.m file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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'))
