THESE/codes/wave/allcode/denoising_dtdwt.m

   1 function y = denoising_dtdwt(x)
   2 % Local Adaptive Image Denoising Algorithm
   3 % Usage :
   4 %        y = denoising_dtdwt(x)
   5 % INPUT :
   6 %        x - a noisy image
   7 % OUTPUT :
   8 %        y - the corresponding denoised image
   9
  10 % Set the windowsize and the corresponding filter
  11 windowsize  = 7;
  12 windowfilt = ones(1,windowsize)/windowsize;
  13
  14 % Number of Stages
  15 J = 6;
  16 I=sqrt(-1);
  17
  18 % symmetric extension
  19 L = length(x); % length of the original image.
  20 N = L+2^J;     % length after extension.
  21 x = symextend(x,2^(J-1));
  22
  23 load nor_dualtree    % run normaliz_coefcalc_dual_tree to generate this mat file.
  24
  25 % Forward dual-tree DWT
  26 % Either FSfarras or AntonB function can be used to compute the stage 1 filters
  27 %[Faf, Fsf] = FSfarras;
  28 [Faf, Fsf] = AntonB;
  29 [af, sf] = dualfilt1;
  30 W = cplxdual2D(x, J, Faf, af);
  31 W = normcoef(W,J,nor);
  32
  33
  34 % Noise variance estimation using robust median estimator..
  35 tmp = W{1}{1}{1}{1};
  36 Nsig = median(abs(tmp(:)))/0.6745;
  37
  38 for scale = 1:J-1
  39     for dir = 1:2
  40         for dir1 = 1:3
  41
  42             % Noisy complex coefficients
  43             %Real part
  44             Y_coef_real = W{scale}{1}{dir}{dir1};
  45             % imaginary part
  46             Y_coef_imag = W{scale}{2}{dir}{dir1};
  47             % The corresponding noisy parent coefficients
  48             %Real part
  49             Y_parent_real = W{scale+1}{1}{dir}{dir1};
  50             % imaginary part
  51             Y_parent_imag = W{scale+1}{2}{dir}{dir1};
  52             % Extend noisy parent matrix to make the matrix size the same as the coefficient matrix.
  53             Y_parent_real  = expand(Y_parent_real);
  54             Y_parent_imag   = expand(Y_parent_imag);
  55
  56             % Signal variance estimation
  57             Wsig = conv2(windowfilt,windowfilt,(Y_coef_real).^2,'same');
  58             Ssig = sqrt(max(Wsig-Nsig.^2,eps));
  59
  60             % Threshold value estimation
  61             T = sqrt(3)*Nsig^2./Ssig;
  62
  63             % Bivariate Shrinkage
  64             Y_coef = Y_coef_real+I*Y_coef_imag;
  65             Y_parent = Y_parent_real + I*Y_parent_imag;
  66             Y_coef = bishrink(Y_coef,Y_parent,T);
  67             W{scale}{1}{dir}{dir1} = real(Y_coef);
  68             W{scale}{2}{dir}{dir1} = imag(Y_coef);
  69
  70         end
  71     end
  72 end
  73
  74 % Inverse Transform
  75 W = unnormcoef(W,J,nor);
  76 y = icplxdual2D(W, J, Fsf, sf);
  77
  78 % Extract the image
  79 ind = 2^(J-1)+1:2^(J-1)+L;
  80 y = y(ind,ind);