% Usage:  Estimate = IHT_Jarvis(y,x,w1,w2,L)
%
% Created by Ramesh Neelamani, Rice University on 11/9/99
%
% Output: Estimate : Inverse halftoned image 
%
% Inputs:
%         y  : half-toned with Jarvis error diffusion
%         w1 : wavelet filter for hard thresholding (reference)
%                                        eg. daubcqf(4,'min');
%         w2 : wavelet filter for  wdwf. eg. daubcqf(6,'min');
%         L  : Number of wavelet levels. eg. 4.
%         x  : original signal (just used to calculate the PSNR)
%
% If the original signal is not provided, then PSNR cannot be calculated.
% If w1,w2,L are not provided, then they are set to default values. 
% 
% This function performs inverse halftoning of halftones obtained with
% Jarvis error diffusion. 
% 
% 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(y,x,w1,w2,L)

M = length(y);
if (nargin == 1 | nargin == 2)
  w1 = daubcqf(4,'min');
  w2 = daubcqf(6,'min');
  L = min(4,log2(M));
end

DisplayVar = 1;
if nargin == 1
  DisplayVar = 0; % Will not calculate PSNR
end

h = zeros(M,M);
K = 4.45;

% FFts of original signal and halftoned signal.
Y = fft2(y);

% Jarvis filter
jar = zeros(M,M);
jar(1,2) = 7;
jar(1,3) = 5;
jar(2,1) = 7;
jar(2,2) = 5;
jar(2,3) = 3;
jar(3,1) = 5;
jar(3,2) = 3;
jar(3,3) = 1;
jar(2,M) = 5;
jar(2,M-1) = 3;
jar(3,M) = 3;
jar(3,M-1) = 1;
h  = jar/48;

H = fft2(h);

% Transfer function for the signal Y = PX + QN  

P = K./(1 + (K-1)*H);
Q = (1-H)./(1+ (K-1)*H);
 
q_by_p = real(ifft2(Q./P));
clear H h Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WInHD begins
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[q_by_pw, q_by_ps] = mrdwt_cycle2D(q_by_p,w1,L);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note that i have re used the same variables to save memory utilization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% getting the multiplicative constants for the noise variance at
% different wavelet levels

for i = 1:L
  for j = 1:3
    mult_const(i,j) = matrix_norm(q_by_pw(:,:,i,j),2);
  end
end

% step 1: inversion
x_noisy = real(ifft2(Y./P));

% Hard thresholding
[x_noisyw, x_noisys] = mrdwt_cycle2D(x_noisy,w1,L);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate sigma of noise N using standard deviation in finest scale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
temp = (x_noisyw(:,:,L,3));
sigma = std(temp(:))/(matrix_norm(q_by_pw(:,:,L,3),2)); % division is to
                                                        % obtain white
                                                        % noise variance.


ThFactor = 2.75;
for i = 1:L
  for j = 1:3
    x_noisyw(:,:,i,j) = HardTh(x_noisyw(:,:,i,j), ThFactor*sigma*mult_const(i,j));
  end
end
x_noisys = HardTh(x_noisys, ThFactor*sigma*matrix_norm(q_by_ps,2));

x_winhd = mirdwt_cycle2D(x_noisyw,x_noisys,w1,L); % this is the intial
                                                % rough estimate
clear x_noisyw X_noisys

% Wavelet domain wiener filtering
[q_by_pw, q_by_ps] = mrdwt_cycle2D(q_by_p,w2,L);

for i = 1:L
  for j = 1:3
    mult_const(i,j) = matrix_norm(q_by_pw(:,:,i,j),2);
  end
end
temp = sigma * matrix_norm(q_by_ps,2);
clear q_by_pw q_by_ps q_by_p

% Wiener filtering
[x_noisyw, x_noisys] = mrdwt_cycle2D(x_noisy,w2,L);
[x_refw, x_refs] = mrdwt_cycle2D(x_winhd,w2,L);

temp_sigma = sigma*mult_const;
for i = 1:L
  for j = 1:3
    x_noisyw(:,:,i,j) = x_noisyw(:,:,i,j).* ...
	( x_refw(:,:,i,j).^2 ./ ((temp_sigma(i,j))^2 + x_refw(:,:,i,j).^2));
  end
end

x_noisys = x_noisys .* ( x_refs.^2 ./ (temp^2 + x_refs.^2));
clear x_refs xrefw
x_wdwf = mirdwt_cycle2D(x_noisyw,x_noisys,w2,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  with Jarvis')

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

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



