dt_bivashrink23.m

基于小波变换的多种去噪方法在matlab上的实现

function y = dt_BivaShrink23(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;

windowfilt2=[1 1 1;1 0 1;1 1 1]/8;

% 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));

% Forward dual-tree DWT
% Either FSfarras or AntonB function can be used to compute the stage 1 filters  
[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 size the same as the coefficient matrix.
            Y_parent_real  = expand(Y_parent_real);
            Y_parent_imag   = expand(Y_parent_imag); 
            
            %----计算邻域系数值----
            Y_adjacent_real=sqrt(conv2((Y_coef_real).^2,windowfilt2,'same'));
            Y_adjacent_imag=sqrt(conv2((Y_coef_imag).^2,windowfilt2,'same'));
                        
            % Signal variance estimation
            Wsig = conv2(windowfilt,windowfilt,(Y_coef_real).^2,'same');%---用Gauss分布的ML估计系数方差----
            Ssig = sqrt(max(Wsig-Nsig.^2,eps));          
            % Threshold value estimation
            T = sqrt(3)*Nsig^2./Ssig;            
            % 
            Y_coef = Y_coef_real+I*Y_coef_imag;
            Y_parent = Y_parent_real + I*Y_parent_imag;
            Y_adjacent = Y_adjacent_real + I*Y_adjacent_imag;            
            
            % Bivariate Shrinkage---邻域系数,父层系数---
            R  = sqrt(abs(Y_parent).^2 + abs(Y_adjacent).^2);
            R = R - T;
            R  = R .* (R > 0);
            Y_coef = Y_coef .* R./(R+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);