THESE/codes/wave/allcode/denoising_dwt.m

   1 function y = denoising_dwt(x)
   2 % Local Adaptive Image Denoising Algorithm
   3 % Usage :
   4 %        y = denoising_dwt(x)
   5 % INPUT :
   6 %        x - a noisy image
   7 % OUTPUT :
   8 %        y - the corresponding denoised image
   9
  10 % Adjust windowsize and the corresponding filter
  11 windowsize  = 7;
  12 windowfilt = ones(1,windowsize)/windowsize;
  13
  14 % Number of Stages
  15 L = 6;
  16
  17 % symmetric extension
  18 N = length(x);
  19 N = N+2^L;
  20 x = symextend(x,2^(L-1));
  21
  22 % forward transform
  23 [af, sf] = farras;
  24 W = dwt2D(x,L,af);
  25
  26 % Noise variance estimation using robust median estimator..
  27 tmp = W{1}{3};
  28 Nsig = median(abs(tmp(:)))/0.6745;
  29
  30 for scale = 1:L-1
  31     for dir = 1:3
  32
  33         % noisy coefficients
  34         Y_coefficient = W{scale}{dir};
  35
  36         % noisy parent
  37         Y_parent = W{scale+1}{dir};
  38
  39         % extent Y_parent to make the matrix size be equal to Y_coefficient
  40         Y_parent = expand(Y_parent);
  41
  42         % Signal variance estimation
  43
  44         Wsig = conv2(windowfilt,windowfilt,(Y_coefficient).^2,'same');
  45         Ssig = sqrt(max(Wsig-Nsig.^2,eps));
  46
  47         % Threshold value estimation
  48         T = sqrt(3)*Nsig^2./Ssig;
  49
  50         % Bivariate Shrinkage
  51         W{scale}{dir} = bishrink(Y_coefficient,Y_parent,T);
  52
  53     end
  54 end
  55
  56
  57 % Inverse Transform
  58 y = idwt2D(W,L,sf);
  59
  60 % Extract the image
  61 y = y(2^(L-1)+1:2^(L-1)+512,2^(L-1)+1:2^(L-1)+512);