BivaShrink123.m
上传用户:sla11nk8
上传日期:2013-03-09
资源大小:21k
文件大小:2k
源码类别:

其他

开发平台:

Matlab

  1. function y = BivaShrink123(x)
  2. % Local Adaptive Image Denoising Algorithm
  3. %---同时用到父层系数和邻域系数信息,是一种三维的萎缩函数---
  4. % Usage :
  5. %        y = BivaShrink123(x)
  6. % INPUT :
  7. %        x - a noisy image
  8. % OUTPUT :
  9. %        y - the corresponding denoised image
  10. % Adjust windowsize and the corresponding filter
  11. windowsize  = 7;
  12. windowfilt = ones(1,windowsize)/windowsize;
  13. windowfilt2=[1 1 1;1 0 1;1 1 1]/8;
  14. % Number of Stages
  15. L = 6;
  16. % symmetric extension
  17. N = length(x);
  18. N = N+2^L;
  19. x = symextend(x,2^(L-1));
  20. % forward transform
  21. [af, sf] = farras;
  22. W = dwt2D(x,L,af); 
  23. % Noise variance estimation using robust median estimator..
  24. tmp = W{1}{3};
  25. Nsig = median(abs(tmp(:)))/0.6745;
  26. for scale = 1:L-1
  27.     for dir = 1:3
  28.         
  29.         % noisy coefficients 
  30.         Y_coefficient = W{scale}{dir};
  31.         
  32.         % noisy parent        
  33.         Y_parent = W{scale+1}{dir};
  34.         
  35.         % extent Y_parent to make the matrix size be equal to Y_coefficient         
  36.         Y_parent = expand(Y_parent);
  37.         
  38.         %----计算邻域系数值----
  39.         Y_adjacent=conv2((Y_coefficient).^2,windowfilt2,'same');
  40.         Y_adjacent=sqrt(Y_adjacent);
  41.         
  42.         % Signal variance estimation
  43.         
  44.         Wsig = conv2(windowfilt,windowfilt,(Y_coefficient).^2,'same');%---用Gauss分布的ML估计系数方差----
  45.         %Wsig_ = conv2(windowfilt,windowfilt,abs(Y_coefficient),'same');%---用Laplace分布的ML估计系数方差----
  46.         %Wsig=(Wsig_*sqrt(2)).^2;
  47.         Ssig = sqrt(max(Wsig-Nsig.^2,eps));
  48.         
  49.         % Threshold value estimation 
  50.         T = sqrt(3)*Nsig^2./Ssig;
  51.         
  52.         % trivariate Shrinkage
  53.         W{scale}{dir} = trishrink1(Y_coefficient,Y_adjacent,Y_parent,T);        
  54.     end
  55. end
  56. % Inverse Transform
  57. y = idwt2D(W,L,sf);
  58. % Extract the image
  59. y = y(2^(L-1)+1:2^(L-1)+512,2^(L-1)+1:2^(L-1)+512);