waveletfilter.m

小波阈值去噪

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% 对图像的小波滤波,这是主程序 ,很多参考资料对小波滤波写的比较简单,即使写的,
%%%%%%%%%一般也只写了一维信号的去噪,并且没有代码,不容易上手.为此我根据自己对小波的理解
%%%%%%%%%对小波边缘去噪的方法,先用matlab写了小波滤波的程序.如果有不对的地方,请指教!

%%%%%%%%%作者:何雄飞
%%%%%%%%%邮箱:[e-mail]

%%%%%%%%%本例图像小波滤波的过程:
%%%%%%%%%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,:));       %%%%%%%%第三层下的噪点个数