2.2 Matlab Implementation of Wavelet-Based Denoising using the Dual-Tree DWT

The implementation of the denoising algorithm is similar to the separable DWT case. There are slight differences since we apply the bivariate shrinkage rule to the magnitudes of the complex coefficients. They are nearly shift invariant, i.e. small signal shifts do not affect the magnitudes, however, they do affect the real and imaginary parts. Therefore magnitude information is more reliable than either real or the imaginary parts. The algorithm does not affect the angle, i.e. we keep it as it is.

The function main_dtdwt.m loads the noisy image, calls the denoising routine and calculates the PSNR value of denoised image.

% Main function
% Usage :
%        main_dtdwt
% INPUT :
%        Raw Barbara image
% OUTPUT :
%        PSNR value of the denoised image
%
% Load clean image
fid = fopen('barbara','r');
s  = fread(fid,[512 512],'unsigned char');
fclose(fid);
N = 512;

% Noise variance
sigma_n = 25;
n = sigma_n*randn(N);

% Add noise
x = s + n;

% Run local adaptive image denoising algorithm using dual-tree DWT.
y = denoising_dtdwt(x);

% Calculate the error
err = s - y;

% Calculate the PSNR value
PSNR = 20*log10(256/std(err(:)))

2. Programs for the Denoising Algorithm

After loading the input image by the "main_dtdwt.m" function, the calculations for the local adaptive image denoising are done by a Matlab function denoising_dtdwt.m. This function calls several subfunctions. The implementation can be summarized as :

1. Set the window size. The signal variance of a coefficent will be estimated using neighboring coefficients in a rectangular region with this window size (see [2] for details).
2. Set how many stages will be used for the wavelet transform.
3. Extend the noisy image. The noisy image will be extended using symetric extension in order to improve the boundary problem with the function symextend.m.
4. Calculate the forward dual-tree DWT using cplxdual2D.m.
5. Estimate the noise variance. The noise variance will be calculated using the robust median estimator.
6. Process each subband seperately in a loop. First the real and imaginary parts of the coefficents and the corresponding parent matrices are prepared for each subband, and the matrices corresponding to the real and imaginary parts of the parent matrix are expanded using a function expand.m in order to make the matrix size the same as the coefficient matrix.
7. Estimate the signal variance and the threshold value : The signal variance for each coefficent is estimated using the window size (see [2] for details), and the threshold value for each coefficient will be calculated and stored in a matrix with the same size as the coefficent matrix.
8. Estimate the magnitude of the complex coefficients. The coefficients will be estimated using the magnitudes of the complex coefficient, its parent and the threshold value with a Matlab function bishrink.m.
9. Calculate the inverse wavelet transform using icplxdual2D.m.
10. Extract the image. The neccessary part of the final image is extracted in order to reverse the symetrical extension.

function y = denoising_dtdwt(x)
% Local Adaptive Image Denoising Algorithm
% Usage :
%        y = denoising_dtdwt(x)
% INPUT :
%        x - a noisy image
% OUTPUT :
%        y - the corresponding denoised image

% Set the windowsize and the corresponding filter
windowsize  = 7;
windowfilt = ones(1,windowsize)/windowsize;

% Number of Stages
J = 6;
I=sqrt(-1);

% symmetric extension
L = length(x); % length of the original image.
N = L+2^J;     % length after extension.
x = symextend(x,2^(J-1));

load nor_dualtree    
% run normaliz_coefcalc_dual_tree to generate this mat file.

% Forward dual-tree DWT
% Either FSfarras or AntonB function can be used for stage 1
%[Faf, Fsf] = FSfarras;
[Faf, Fsf] = AntonB;
[af, sf] = dualfilt1;
W = cplxdual2D(x, J, Faf, af);
W = normcoef(W,J,nor);


% Noise variance estimation using robust median estimator..
tmp = W{1}{1}{1}{1};
Nsig = median(abs(tmp(:)))/0.6745;

for scale = 1:J-1
    for dir = 1:2
        for dir1 = 1:3

            % Noisy complex coefficients
            %Real part
            Y_coef_real = W{scale}{1}{dir}{dir1};
            % imaginary part
            Y_coef_imag = W{scale}{2}{dir}{dir1};
            % The corresponding noisy parent coefficients
            %Real part
            Y_parent_real = W{scale+1}{1}{dir}{dir1};
            % imaginary part
            Y_parent_imag = W{scale+1}{2}{dir}{dir1};
            % Extend noisy parent matrix to make the matrix the
            % same size as the coefficient matrix.
            Y_parent_real  = expand(Y_parent_real);
            Y_parent_imag   = expand(Y_parent_imag);

            % Signal variance estimation
            Wsig = conv2(windowfilt,windowfilt,(Y_coef_real).^2,'same');
            Ssig = sqrt(max(Wsig-Nsig.^2,eps));

            % Threshold value estimation
            T = sqrt(3)*Nsig^2./Ssig;

            % Bivariate Shrinkage
            Y_coef = Y_coef_real+I*Y_coef_imag;
            Y_parent = Y_parent_real + I*Y_parent_imag;
            Y_coef = bishrink(Y_coef,Y_parent,T);
            W{scale}{1}{dir}{dir1} = real(Y_coef);
            W{scale}{2}{dir}{dir1} = imag(Y_coef);

        end
    end
end

% Inverse Transform
W = unnormcoef(W,J,nor);
y = icplxdual2D(W, J, Fsf, sf);

% Extract the image
ind = 2^(J-1)+1:2^(J-1)+L;
y = y(ind,ind);
 

3. Summary of programs:

1. AntonB.m - Compute filters for stage 1
2. dualfilt1.m - Compute filters for stage the remaining stages
3. bishrink.m - Estimation using the bivariate shrinkage function
4. expand.m - Expand the parent matrices
5. cplxdual2D.m - Forward dual-tree DWT
6. icplxdual2D.m - Inverse dual-tree DWT
7. symextend.m - Extend the image symmetrically
8. denoising_dtdwt.m - Local adaptive image processing routine
9. main_dtdwt.m - Main function