% Denoising by UWT, undecimated wavelet transform

  % Read image;
img = double(imread('einstein.bmp'));
[N,J] = quadlength(img);
% Add white Gaussian noise.
imn = img + 20*randn(size(img));

% Coarsest decomposition scale.
coarsest = 2;

% Generate Symmlet 4 CMF Filter.
qmf = MakeONFilter('Symmlet',4);

% Compute UWT and DWT of noisy image.
uwc = FWT2_TI(imn,coarsest,qmf);
wc  = FWT2_PO(imn,coarsest,qmf);

% Estimate noise std from diagonal band at finest scale.
hh = uwc(2*N+1:3*N,:);
sigma = MAD(hh(:));

% Hard-threshold UWT coeffs: finest scale at 4*sigma and other scales at 3*sigma.
uwc(1:3*N,:)       = HardThresh(uwc(1:3*N,:),4*sigma);       % Finest.
uwc(3*N+1:end-N,:) = HardThresh(uwc(3*N+1:end-N,:),3*sigma); % Other scales.

% Hard-threshold DWT coeffs:
ll  = wc(1:coarsest,1:coarsest);
wc(N/2+1:N,1:N/2)   = HardThresh(wc(N/2+1:N,1:N/2),4*sigma);   % Finest vertical.
wc(1:N/2,N/2+1:N)   = HardThresh(wc(1:N/2,N/2+1:N),4*sigma);   % Finest horizontal.
wc(N/2+1:N,N/2+1:N) = HardThresh(wc(N/2+1:N,N/2+1:N),4*sigma); % Finest diagonal.
wc(1:N/2,1:N/2)     = HardThresh(wc(1:N/2,1:N/2),3*sigma);     % Other scales.
wc(1:coarsest,1:coarsest) = ll;

% Reconstruct.
imdenuwt = IWT2_TI(uwc,coarsest,qmf);
imdendwt = IWT2_PO(wc,coarsest,qmf);

% Display.
subplot(221)
     imagesc(img);axis image;axis off
set(gca,'FontSize',14);
title('(a)');
subplot(222)
     imagesc(imn);axis image;axis off
set(gca,'FontSize',14);
title('(b)');
subplot(223)
     imagesc(imdenuwt);axis image;axis off
set(gca,'FontSize',14);
title('(c)');
subplot(224)
     imagesc(imdendwt);axis image;axis off
set(gca,'FontSize',14);
title('(d)');
colormap('gray')