% This function performs 2-D wavelet denoising of noisy images using  
% Hyperbolic shrinkage with Donoho's universal" threshold.
% Inputs:  s -- original grey-scale image normalized to [0, 1] 
%          x -- noise-corrupted image, obtained as
%               randn('state',7); 
%               x = s + 0.1 * randn(size(s));
%         h0 -- impulse response of the lowpass analysis filter
%         h1 -- impulse response of the highpass analysis filter
%         f0 -- impulse response of the lowpass synthesis filter
%         f1 -- impulse response of the highpass synthesis filter
% Output: xh -- denoised image (normalized)
% Written by W.-S. Lu, University of Victoria.
% Last modified: March 3, 2003.

function xh = deno_2(s,x,h0,h1,f0,f1)

sz = size(x);
N = min(sz);
p = ceil(log10(N)/log10(2));

% 2-D Wavelet Decomposition

y = wt2d(x,h0,h1,p);

% Estimate Noise Variance

d1 = y(N/2+1:N,N/2+1:N);
db = mean(mean(d1));
s2 = sum(sum((d1-db).^2))/((N/2)^2-1);

% Compute Donoho's Universal Threshold

delta = 2*sqrt(log10(N)*s2);

% Hyperbolic Shrinkage of Wavelet Coefficients xd

yh = sign(y).*sqrt(max(y.^2-delta^2,zeros(N,N)));

% Reconstruct the Denoised Image

xh = iwt2d(yh,f0,f1,p);
 
% Project xh onto [0, 1]

xh = min(ones(size(xh)),xh);
xh = max(zeros(size(xh)),xh);

% Performancve Evaluation

nf = norm(s,'fro');

% Signal-to-Noise ratio Before Denoising:

SNR1 = 20*log10(nf/norm(x-s,'fro'));

% Signal-to-Noise Ratio After Denoising:

SNR2 = 20*log10(nf/norm(xh-s,'fro'));

% Signal-to-Noise Improvement

ISNR = SNR2-SNR1;
disp('before-after-improvement SNR:')
[SNR1 SNR2 ISNR]

% Display Results

map = gray(256);
figure(1)
subplot(221)
imshow(s*255,map)
title('(a) origimal image')
axis('square')
subplot(222)
imshow(x*255,map)
title('(b) noise-corrupted image')
axis('square')
subplot(223)
imshow(xh*255,map)
title('(c) denoised image')
axis('square')