2.1 Matlab Implementation of Wavelet-Based Denoising using the Separable DWT

In our implementation, the main function calls the algorithm as a function. This function loads the noisy image, calls the denoising routine and calculates the PSNR value of the denoised image. This main function is implemented with the Matlab function main_dwt.m

% Main function
% Usage :
%        main_dwt
% INPUT :
%        Raw Lena image
% OUTPUT :
%        PSNR value of the denoised image
%
% Load clean image
fid = fopen('lena','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
y = denoising_dwt(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_dwt.m" function, the calculations for the local adaptive image denoising are done by a Matlab function denoising_dwt.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 reduce the boundary problem with a function symextend.m.
4. Calculate the forward wavelet transform.
5. Estimate the noise variance. The noise variance will be calculated using the robust median estimator.
6. Process each subband seperately. We will process each subband in a loop since our implementation stores the wavelet coefficients in a cell array. First the coefficent and the corresponding parent matrices are prepared for each subband, and the parent matrix is expanded using the 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 coefficients. The coefficients will be estimated using the noisy coefficient, its parent, and the estimated threshold value with the Matlab function bishrink.m.
9. Calculate the inverse wavelet transform.
10. Extract the image. The neccessary part of the final image is extracted in order to reverse the symetric extension operation.

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

% Adjust windowsize and the corresponding filter
windowsize  = 7;
h = ones(1,windowsize)/windowsize;

% Number of Stages
L = 6;

% symmetric extension
N = length(x);
N = N+2^L;
x = symextend(x,2^(L-1));

% forward transform
[af, sf] = farras;
W = dwt2D(x,L,af); 

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

for scale = 1:L-1
    for dir = 1:3
        
        % noisy coefficients 
        Y_coefficient = W{scale}{dir};
        
        % noisy parent        
        Y_parent = W{scale+1}{dir};
        
        % extend Y_parent to make it the same size as Y_coefficient
        Y_parent = expand(Y_parent);
        
        % Signal variance estimation
        
        Wsig = conv2(h,h,(Y_coefficient).^2,'same');
        Ssig = sqrt(max(Wsig-Nsig.^2,eps));
        
        % Threshold value estimation 
        T = sqrt(3)*Nsig^2./Ssig;
        
        % Bivariate Shrinkage
        W{scale}{dir} = bishrink(Y_coefficient,Y_parent,T);
        
    end
end


% Inverse Transform
y = idwt2D(W,L,sf);

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

3. Summary of programs:

1. farras.m - A set of perfect reconstruction filters
2. bishrink.m - Estimation using the bivariate shrinkage function
3. expand.m - Expand the parent matrices
4. dwt.m - Forward wavelet transform
5. idwt.m - Inverse wavelet transform
6. symextend.m - Extend the image symmetrically
7. denoising_dwt.m - Local adaptive image processing routine
8. main_dwt.m - Main_dwt.m