%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% 對圖像的小波濾波,這是主程序 ,很多參考資料對小波濾波寫的比較簡單,即使寫的,
%%%%%%%%%一般也只寫了一維信號的去噪,并且沒有代碼,不容易上手.為此我根據自己對小波的理解
%%%%%%%%%對小波邊緣去噪的方法,先用matlab寫了小波濾波的程序.如果有不對的地方,請指教!

%%%%%%%%%作者:何雄飛
%%%%%%%%%郵箱:steghexiongfei@163.com

%%%%%%%%%本例圖像小波濾波的過程:
%%%%%%%%%1)讀取圖像,空域,這里采用亮度系數
%%%%%%%%%2)對數據進行小波變換,三層haar小波變換,因為haar小波比較簡單,也可以用JPEG那個小波變換
%%%%%%%%%3)估計噪聲系數,分三個方向,水平H,垂直V,對角D.
%%%%%%%%%4)從第三層開始小波變換,權重系數weight和threshold可以參考部分資料,因為噪聲對三層的影響不一樣
%%%%%%%%%5)只有在大尺寸是邊緣的在小尺寸的時候才是邊緣
%%%%%%%%%6)從3->1開始重構圖像
%%%%%%%%%7)可以看到噪聲被提取出來.


function  WA0=wavefilter(imfile)

      Image = imread(imfile); 
      [r,l,color]=size(Image);
      if color==3
          Image=Image(:,:,1);
      end
      A0=double(Image);
      %%%%%%%%%%%%% Three-level haar wavelet %%%%%%%%%%%%%%%%%%%%%%%%
      [A1,H1,V1,D1] = dwt2(double(A0),'Haar');
      [A2,H2,V2,D2] = dwt2(double(A1),'Haar');
      [A3,H3,V3,D3] = dwt2(double(A2),'Haar');
      %%%%%%%%%%%%初始化小波系數系數  備份,后面要用到!
      WA0=A0;
      
      WA1=A1;WH1=H1;WV1=V1;WD1=D1;
      [r1,l1]=size(WA1);
      %%%%%%%%%%%%%%%
      WA2=A2;WH2=H2;WV2=V2;WD2=D2;
      [r2,l2]=size(WA2);
      %%%%%%%%%%%%%%%%
      WA3=A3;WH3=H3;WV3=V3;WD3=D3;
      [r3,l3]=size(WA3);
      %%%%%%%%%%%%%%%%%%%初始化掩膜，即屏蔽濾波器MASK=0
      MaskH1=zeros(r1,l1); MaskV1=zeros(r1,l1); MaskD1=zeros(r1,l1);
      
      MaskH2=zeros(r2,l2); MaskV2=zeros(r2,l2); MaskD2=zeros(r2,l2);
      
      MaskH3=zeros(r3,l3); MaskV3=zeros(r3,l3); MaskD3=zeros(r3,l3);
      %%%%%%%%%%%%%%%%%%引入與尺度相關的權重系數   
      weight=zeros(1,3);
      weight(1)=1.15; weight(2)=1.06; weight(3)=1.00;
      
      %%%%%%%%%%%%%%%%引入噪聲能量的權值   可以適當修改
      threshold=zeros(3,3);
      threshold(1,1)=1.10;threshold(1,2)=1.10;threshold(1,3)=1.10;
      threshold(2,1)=1.30;threshold(2,2)=1.30;threshold(2,3)=1.30;
      threshold(3,1)=1.50;threshold(3,2)=1.50;threshold(3,3)=1.50;
      
      %%%%%%%%%%%%%%%%初始化相關系數
      CorH1=zeros(r1,l1); CorV1=zeros(r1,l1); CorD1=zeros(r1,l1);
      
      CorH2=zeros(r2,l2); CorV2=zeros(r2,l2); CorD2=zeros(r2,l2);
      
      CorH3=zeros(r3,l3); CorV3=zeros(r3,l3); CorD3=zeros(r3,l3);
      %%%%%%%%%%%%%%計算不同尺度下小波系數的相關性%%%%%%%%%%%%%
      CorH1=WH1.*WA1;CorV1=WV1.*WA1;CorD1=WD1.*WA1;   %%計算細節和主體之間的相關性，分3個子帶
     
      CorH2=WH2.*WA2;CorV2=WV2.*WA2;CorD2=WD2.*WA2;
       
      CorH3=WH3.*WA3;CorV3=WV3.*WA3;CorD3=WD3.*WA3;
      
      %%%%%%%%%%%%%%%初始化邊緣點個數
      N1=r1*l1; N2=r2*l2; N3=r3*l3;   %%%%%%%%%%%%%%總的點的個數
      k=zeros(3,3);                   %%%%%%%%%%%%%邊緣點的個數
      
      %%%%%%%%%尋找掩膜的寫成一個函數來調用  此處先求     H  
      pw=sum(sum(WH1.*WH1));
      pcor=sum(sum(CorH1.*CorH1));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r1,l1);
      NCor=CorH1*normalization;
      for i=1:1:r1
          for j=1:1:l1
              if abs(NCor(i,j))>=abs(WH1(i,j))
                  CorH1(i,j)=0;
                  WH1(i,j)=0;
                  MaskH1(i,j)=1;
                  k(1,1)=k(1,1)+1;
              end
          end
      end
      %%%%%%%%%%%%利用一次小波變換的一次比較，剩余的小波系數作為噪聲
      %%%%%%%%%%%%然后，利用這個噪聲推廣到全局。
      pw=sum(sum(WH1.*WH1));
      SNH=pw;
      snh1=SNH/(N1-k(1,1));
      snh2=snh1/4;
      snh3=snh1/16;
      %%%%%%%%%%開始倒著計算掩膜，這樣是改進的方法！    H3
      pw=sum(sum(WH3.*WH3));
      pcor=sum(sum(CorH3.*CorH3));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r3,l3);
      NCor=CorH3*normalization;     
      while (0.95*pw>threshold(3,1)*(N3-k(3,1))*snh1)&(k(3,1)<N3) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(3,1);
         [WH3,CorH3,MaskH3,k(3,1),pw,NCor]= findmask(WH3,CorH3,MaskH3,k(3,1),NCor,weight(3));
         if k(3,1)==a
            break;
         end
      end
      %%%%%%%%%%%%%開始計算第3層水平小波變換           H2
      pw=sum(sum(WH2.*WH2));
      pcor=sum(sum(CorH2.*CorH2));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r2,l2);
      NCor=CorH2*normalization;     
      while (0.95*pw>threshold(2,1)*(N2-k(2,1))*snh1)&(k(2,1)<N2) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(2,1);
         [WH2,CorH2,MaskH2,k(2,1),pw,NCor]= findmask(WH2,CorH2,MaskH2,k(2,1),NCor,weight(2));
         if k(2,1)==a
            break;
         end
      end
      %%%%%%%%%%%%%%%改進的算法，也需要計算第一層，因為引入了權重和百分比，還要由此計算噪聲掩膜
      pw=sum(sum(WH1.*WH1));
      pcor=sum(sum(CorH1.*CorH1));                  % H1
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r1,l1);
      NCor=CorH1*normalization; 
      while (0.95*pw>threshold(1,1)*(N1-k(1,1))*snh1)&(k(1,1)<N1) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(1,1);
         [WH1,CorH1,MaskH1,k(1,1),pw,NCor]= findmask(WH1,CorH1,MaskH1,k(1,1),NCor,weight(1));
         if k(1,1)==a
            break;
         end
      end
      %%%%%%%%%%即小波濾波的第三個要素，只有是大尺度的邊緣，才能是小尺度的邊緣。
      MaskH2=MaskH2.*pixeldup(MaskH3,1);
      MaskH1=MaskH1.*pixeldup(MaskH2,1);
      
      %%%%%%%%%尋找掩膜的寫成一個函數來調用  此處先求          V  
      pw=sum(sum(WV1.*WV1));
      pcor=sum(sum(CorV1.*CorV1));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r1,l1);
      NCor=CorV1*normalization;
      for i=1:1:r1
          for j=1:1:l1
              if abs(NCor(i,j))>=abs(WV1(i,j))
                  CorV1(i,j)=0;
                  WV1(i,j)=0;
                  MaskV1(i,j)=1;
                  k(1,2)=k(1,2)+1;
              end
          end
      end
      %%%%%%%%%%%%利用一次小波變換的一次比較，剩余的小波系數作為噪聲
      %%%%%%%%%%%%然后，利用這個噪聲推廣到全局。
      pw=sum(sum(WV1.*WV1));
      SNV=pw;
      snv1=SNV/(N1-k(1,2));
      snv2=snv1/4;
      snv3=snv1/16;
      %%%%%%%%%%開始倒著計算掩膜，這樣是改進的方法！    V3
      pw=sum(sum(WV3.*WV3));
      pcor=sum(sum(CorV3.*CorV3));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r3,l3);
      NCor=CorV3*normalization;     
      while (0.95*pw>threshold(3,2)*(N3-k(3,2))*snv1)&(k(3,2)<N3) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(3,2);
         [WV3,CorV3,MaskV3,k(3,2),pw,NCor]= findmask(WV3,CorV3,MaskV3,k(3,2),NCor,weight(3));
         if k(3,2)==a
            break;
         end
      end
      %%%%%%%%%%%%%開始計算第3層水平小波變換           V2
      pw=sum(sum(WV2.*WV2));
      pcor=sum(sum(CorV2.*CorV2));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r2,l2);
      NCor=CorV2*normalization;     
      while (0.95*pw>threshold(2,2)*(N2-k(2,2))*snv1)&(k(2,2)<N2) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(2,2);
         [WV2,CorV2,MaskV2,k(2,2),pw,NCor]= findmask(WV2,CorV2,MaskV2,k(2,2),NCor,weight(2));
         if k(2,2)==a
            break;
         end
      end
      %%%%%%%%%%%%%%%改進的算法，也需要計算第一層，因為引入了權重和百分比，還要由此計算噪聲掩膜
      pw=sum(sum(WV1.*WV1));
      pcor=sum(sum(CorV1.*CorV1));                  % V1
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r1,l1);
      NCor=CorV1*normalization; 
      while (0.95*pw>threshold(1,2)*(N1-k(1,2))*snv1)&(k(1,2)<N1) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(1,2);
         [WV1,CorV1,MaskV1,k(1,2),pw,NCor]= findmask(WV1,CorV1,MaskV1,k(1,2),NCor,weight(1));
         if k(1,2)==a
            break;
         end
      end
      %%%%%%%%%%即小波濾波的第三個要素，只有是大尺度的邊緣，才能是小尺度的邊緣。
      MaskV2=MaskV2.*pixeldup(MaskV3,1);
      MaskV1=MaskV1.*pixeldup(MaskV2,1);
      
      %%%%%%%%%尋找掩膜的寫成一個函數來調用  此處先求      D  
      pw=sum(sum(WD1.*WD1));
      pcor=sum(sum(CorD1.*CorD1));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r1,l1);
      NCor=CorD1*normalization;
      for i=1:1:r1
          for j=1:1:l1
              if abs(NCor(i,j))>=abs(WD1(i,j))
                  CorD1(i,j)=0;
                  WD1(i,j)=0;
                  MaskD1(i,j)=1;
                  k(1,3)=k(1,3)+1;
              end
          end
      end
      %%%%%%%%%%%%利用一次小波變換的一次比較，剩余的小波系數作為噪聲
      %%%%%%%%%%%%然后，利用這個噪聲推廣到全局。
      pw=sum(sum(WD1.*WD1));
      SND=pw;
      snd1=SND/(N1-k(1,3));
      snd2=snd1/4;
      snd3=snd1/16;
      %%%%%%%%%%開始倒著計算掩膜，這樣是改進的方法！    D3
      pw=sum(sum(WD3.*WD3));
      pcor=sum(sum(CorD3.*CorD3));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r3,l3);
      NCor=CorD3*normalization;
      while (0.95*pw>threshold(3,3)*(N3-k(3,3))*snd1)&(k(3,3)<N3) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(3,3);
         [WD3,CorD3,MaskD3,k(3,3),pw,NCor]= findmask(WD3,CorD3,MaskD3,k(3,3),NCor,weight(3));
         if k(3,3)==a
            break;
         end
      end
      %%%%%%%%%%%%%開始計算第3層水平小波變換           D2
      pw=sum(sum(WD2.*WD2));
      pcor=sum(sum(CorD2.*CorD2));
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r2,l2);
      NCor=CorD2*normalization;     
      while (0.95*pw>threshold(2,3)*(N2-k(2,3))*snd1)&(k(2,3)<N2) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(2,3);
         [WD2,CorD2,MaskD2,k(2,3),pw,NCor]= findmask(WD2,CorD2,MaskD2,k(2,3),NCor,weight(2));
         if k(2,3)==a
            break;
         end
      end
      %%%%%%%%%%%%%%%改進的算法，也需要計算第一層，因為引入了權重和百分比，還要由此計算噪聲掩膜
      pw=sum(sum(WD1.*WD1));
      pcor=sum(sum(CorD1.*CorD1));                  % D1
      normalization=sqrt(pw/(pcor+eps));
      NCor=zeros(r1,l1);
      NCor=CorD1*normalization;  
      while (0.95*pw>threshold(1,3)*(N1-k(1,3))*snd1)&(k(1,3)<N1) %新的約束條件，增加了噪聲影響權重，還有增加了數值百分比
          a=k(1,3);
         [WD1,CorD1,MaskD1,k(1,3),pw,NCor]= findmask(WD1,CorD1,MaskD1,k(1,3),NCor,weight(1));
         if k(1,3)==a
            break;
         end
      end
      %%%%%%%%%%即小波濾波的第三個要素，只有是大尺度的邊緣，才能是小尺度的邊緣。
      MaskD2=MaskD2.*pixeldup(MaskD3,1);
      MaskD1=MaskD1.*pixeldup(MaskD2,1);
            
      %%%%%%%%%%%%%%%%%對原始數據施加屏蔽濾波
      WH1=H1.*MaskH1; WV1=V1.*MaskV1; WD1=D1.*MaskD1;
      
      WH2=H2.*MaskH2; WV2=V2.*MaskV2; WD2=D2.*MaskD2;
      
      WH3=H3.*MaskH3; WV3=V3.*MaskV3; WD3=D3.*MaskD3;
      %%%%%%%%%%%%%%%%重構小波，使得能夠重新恢復圖像
      WA2=idwt2(WA3,WH3,WV3,WD3,'haar');
      WA1=idwt2(WA2,WH2,WV2,WD2,'haar');
      WA0=idwt2(WA1,WH1,WV1,WD1,'haar');   %%%%%%%%重構出濾波后的圖像
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      D=A0-WA0;               %%%%%%%得到噪聲圖像
      N=D.*D;                 %%%%%%%得到噪聲能量
      SUMN=sum(sum(N));       %%%%%%%噪聲總能量
      SUMNOISE=SNH+SNV+SND;
      imshow(uint8(abs(D)));  %%%%%%%%顯示噪聲圖像
      %%%%%%%%%%%%求出噪點個數
      k(1,:)=r1*l1-k(1,:);
      k(2,:)=r2*l2-k(2,:);
      k(3,:)=r3*l3-k(3,:);
      %%%%%%%%%%%%觀察噪聲層噪點能量
      n=zeros(1,3);
      n(1)=sum(k(1,:));       %%%%%%%%小尺度下的噪點個數
      n(2)=sum(k(2,:));       %%%%%%%%第二層下的噪點個數
      n(3)=sum(k(3,:));       %%%%%%%%第三層下的噪點個數