% Usage:  Estimate = IHT_Jarvis_cwt(y,x,L)
%
% Created by Ramesh Neelamani, Rice University on 10/02/02
%
% Output: Estimate : Inverse halftoned image 
%
% Inputs:
%         y  : half-toned with  Jarvis error diffusion
%         L  : Number of wavelet levels
%         x  : original signal (just used to calculate the PSNR)
%
% If the original signal is not provided, then PSNR cannot be calculated.
% 
% This function performs inverse halftoning of halftones obtained with
% Jarvis error diffusion.  This program uses Complex Wavelets.
% 
% 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/
%
 
function [x_wdwf] = IHT_Jarvis_cwt(y,x,L)

DisplayVar = 1;
M = length(y);

if nargin == 1
  DisplayVar = 0; % Will not calculate PSNR
  L = min(4,log2(M));
elseif nargin == 2
   L = min(4,log2(M)); 
end
K = 4.45;

% Floyd-Stienberg filter
jar(1,4) = 7;
jar(1,5) = 5;
jar(2,3) = 7;
jar(2,4) = 5;
jar(2,5) = 3;
jar(3,3) = 5;
jar(3,4) = 3;
jar(3,5) = 1;
jar(2,2) = 5;
jar(2,1) = 3;
jar(3,2) = 3;
jar(3,1) = 1;
h  = jar/48;


% Transfer function for the signal Y = PX + QN  P -> fft2 of p
p_inv = (K-1)*h;
p_inv(1,3) = 1;
p_inv = p_inv/K;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WInHD begins
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% getting the multiplicative constants for the noise variance at
% different wavelet levels
noise_runs = 1;
[varmat1,varmat2] = colored_variances_b(noise_runs,M,L,h,K);

% step 1: inversion of P by convolving with p_inv
temp_y = [y(M,M-3:M) y(M,:);  y(:,M-3:M) y];
x_noisy = conv2(p_inv,temp_y);
x_noisy = x_noisy(2:M+1,7:M+6);

% Hard thresholding
[x_noisyw1, x_noisyw2] = cdwt2_b(x_noisy,L);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate sigma of noise N using standard deviation in finest scale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
temp1 = x_noisyw1(M/2+1:M,M/2+1:M);
temp2 = x_noisyw2(M/2+1:M,M/2+1:M);

sigma = sqrt((var(temp1(:))/((varmat1(M,M))) + var(temp2(:))/((varmat2(M,M))))/2);   

                          % division is to obtain white noise variance.
                          % averaging for more robust estimate

ThFactor = 2.7;
x_noisyw1 = HardTh(x_noisyw1, ThFactor*sigma*sqrt(varmat1));
x_noisyw2 = HardTh(x_noisyw2, ThFactor*sigma*sqrt(varmat2));

x_winhd = icdwt2_b(x_noisyw1,x_noisyw2,L); % this is the intial
                                         % rough estimate
clear x_noisyw1 x_noisyw2

% Wavelet domain wiener filtering
[varmat1,varmat2] = colored_variances(noise_runs,M,L,h,K);
[x_noisyw1, x_noisyw2] = cdwt2(x_noisy,L);
[x_refw1, x_refw2] = cdwt2(x_winhd,L);

x_noisyw1 = x_noisyw1.*(abs(x_refw1).^2 ./ (abs(x_refw1).^2+sigma^2*varmat1));
x_noisyw2 = x_noisyw2.*(abs(x_refw2).^2 ./ (abs(x_refw2).^2+sigma^2*varmat2));

clear x_refw1 xrefw2
x_wdwf = icdwt2(x_noisyw1,x_noisyw2,L);


if DisplayVar == 1
  figure(200)
  clf
  subplot(1,3,1)
  imagesc(x)
  im_range = get(gca,'CLim');
  imset  
  title('Original')

  subplot(1,3,2)
  imagesc(y)
  imset
  title('Halftoned')

  subplot(1,3,3)
  imagesc(x_wdwf,im_range)
  rmse_wdwf = norm(x(:)-x_wdwf(:));
  error_db_wdwf = 20*log10(norm(x(:))/rmse_wdwf);
  psnr_db_wdwf = 20*log10(M*255/rmse_wdwf);
  imset
  title(strcat('Inverse Halftoned,','PSNR = ',num2str(psnr_db_wdwf),'dB'))
else 
  figure(200)
  clf
  subplot(1,2,1)
  imagesc(y)
  imset
  title('Halftoned')

  subplot(1,2,2)
  imagesc(x_wdwf)
  imset
  title('Inverse Halftoned')
end



