2.2 Matlab Implementation of Wavelet-Based Denoising using the Dual-Tree DWTThe 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.
Table 2.1 Matlab function main_dtdwt.m
% 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 :
Table 2.2 Matlab function denoising_dtdwt.m
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:
Separable DWT
1-D DWT2-D DWT 3-D DWT Dual-tree DWT
1-D DWT2-D DWT 3-D DWT Applications
Denoising IDenoising II |