基于纹理模型的空间自适应小波滤波实现

 

 主程序

clear
[X,map]=imread('teacher4-stand2.png');              
  
X=rgb2gray(X);       
%set(gcf,'Position',[0,0,512,512])
%set(gca,'position',[0 0 1 1])
subplot(2,2,1);
imshow(X); title('原始图像');                    
  
% 生成含噪图像并图示  imnoise(I, 'salt & pepper');
  
init=2055615866;         
  
randn('seed',init);        
  
X=double(X);  
  
% 添加随机噪声  
  XX=X;
 XX=X+20*randn(size(X));    
subplot(2,2,2);
imshow(uint8(XX));                
title(' 含噪图像 ');  
XX=trancetion(XX);
[LL,HL,LH,HH]=mydwt2(XX);%第一层二维小波分解
[m,n]=size(LL);
LL=trancetion(LL);
[cA2,cV2,cH2,cD2]=mydwt2(LL);%第二层小波分解
LH1=ghgh(cH2);
HL1=ghgh(cV2);
HH1=ghgh(cD2);
LL1=ghgh(cA2);
LLt=myidwt2(LL1,HL1,LH1,HH1);
LLt=imresize(LLt,[m,n]);
HL=ghgh(HL);
LH=ghgh(LH);
HH=ghgh(HH);
LL=ghgh(LL);
% LL：低频部分分解系数； HL：垂直方向分解系数；对应V
% LH：水平方向分解系数对应H； HH：对角线方向分解系数。
new_img=myidwt2(LLt,HL,LH,HH);
%nen=myidwt2(,HL,LH,HH);
subplot(224);imshow(uint8(new_img));title('(c)二层目标算法去噪图像');
subplot(223);imshow(uint8(LLt));title('(d)第二层层目标算法图像');

function T=ghgh(x1)
% % %目标算法1，以纹理值分为两部分
a=3;cH11=x1;
[col,woe]=size(cH11);
cH11=MirrorImage(double(cH11),a);
new1=ones(col+2*a,woe+2*a);
for i=a+1:col+a
    for j=a+1:woe+a
        med=(2*(cH11(i-1,j)+cH11(i,j-1)+cH11(i,j+1)+cH11(i+1,j))+(cH11(i-1,j-1)+cH11(i-1,j+1)+cH11(i+1,j-1)+cH11(i+1,j+1)))/12;
        new1(i,j)=med;
    end
end
cH11=cH11((a+1):(col+a),(a+1):(woe+a));%
new1=new1((a+1):(col+a),(a+1):(woe+a));%z(x,y)
med1=mean(new1(:));%取中值作为分界
new2=zeros(col,woe);ne1=[];
new3=zeros(col,woe);ne2=[];%两个空矩阵用于存储部分小波系数
w1=1;
w2=1;
for i=1:col
    for j=1:woe
        m=new1(i,j);%以Z（x,y）作为纹理值
        if m<med1
            new2(i,j)=cH11(i,j);
            ne1(w1)=new2(i,j);w1=w1+1;
        else
            new3(i,j)=cH11(i,j);
             ne2(w2)=new3(i,j);w2=w2+1;
        end
    end
end

    % cH1=new2(:);
     va=var(ne1);dia=median(abs(ne1))/0.6745;
    diax2=va-dia^2;
    T1=sqrt(2)*dia^2/sqrt(diax2);
for i=1:col
    for j=1:woe
        cH1=new2(i,j);
       if abs(cH1)<T1
          cH1=0;
       else
           if cH1>=0
              cH1=abs(cH1)-T1;
           else
              cH1=-1*(abs(cH1)-T1); 
           end
       end
     new2(i,j)=cH1;
    end
end
va=var(ne2);dia=median(abs(ne2))/0.6745;
diax2=va-dia^2;
T2=sqrt(2)*dia^2/sqrt(diax2);
for i=1:col
    for j=1:woe
        cH2=new3(i,j);
        if abs(cH2)<T2
           cH2=0;
        else
            if cH2>=0
              cH2=abs(cH2)-T2;
            else
              cH2=-1*(abs(cH2)-T2); 
            end
        end
        new3(i,j)=cH2;
    end
end
T=new2+new3;%两个矩阵相加，返回系数矩阵。

以下借鉴他人的

function [cA,cD] = mydwt(x,lpd,hpd,dim)
% 函数 [cA,cD]=MYDWT(X,LPD,HPD,DIM) 对输入序列x进行一维离散小波分解，输出分解序列[cA,cD]
% 输入参数：x——输入序列；
% lpd——低通滤波器；
% hpd——高通滤波器；
% dim——小波分解级数。
% 输出参数：cA——平均部分的小波分解系数；
% cD——细节部分的小波分解系数。
 
cA=x; % 初始化cA，cD
cD=[];
for i=1:dim
cvl=conv(cA,lpd); % 低通滤波，为了提高运行速度，调用MATLAB提供的卷积函数conv()
dnl=downspl(cvl); % 通过下抽样求出平均部分的分解系数
cvh=conv(cA,hpd); % 高通滤波
dnh=downspl(cvh); % 通过下抽样求出本层分解后的细节部分系数
cA=dnl; % 下抽样后的平均部分系数进入下一层分解
cD=[cD,dnh]; % 将本层分解所得的细节部分系数存入序列cD
end
 
function y=downspl(x);
% 函数 Y=DOWMSPL(X) 对输入序列进行下抽样，输出序列 Y。
% 下抽样是对输入序列取其偶数位，舍弃奇数位。例如 x=[x1,x2,x3,x4,x5]，则 y=[x2,x4].
 
N=length(x); % 读取输入序列长度
M=floor(N/2); % 输出序列的长度是输入序列长度的一半（带小数时取整数部分）
i=1:M;
y(i)=x(2*i);
 

function [LL,HL,LH,HH]=mydwt2(x)
% 函数 MYDWT2() 对输入的r*c维矩阵 x 进行二维小波分解，输出四个分解系数子矩阵[LL,HL,LH,HH]
% 输入参数：x —— 输入矩阵，为r*c维矩阵。https://blog.csdn.net/baidu_35561918/article/details/52421653
% 输出参数：LL,HL,LH,HH —— 是分解系数矩阵的四个相等大小的子矩阵,大小均为 r/2 * c/2 维
% LL：低频部分分解系数； HL：垂直方向分解系数；
% LH：水平方向分解系数； HH：对角线方向分解系数。
 
lpd=[1/2 1/2];hpd=[-1/2 1/2]; % 默认的低通、高通滤波器
[row,col]=size(x); % 读取输入矩阵的大小
 
for j=1:row % 首先对输入矩阵的每一行序列进行一维离散小波分解
     tmp1=x(j,:);
     [ca1,cd1]=mydwt(tmp1,lpd,hpd,1);
     x(j,:)=[ca1,cd1]; % 将分解系数序列再存入矩阵x中，得到[L|H]
end
for k=1:col % 再对输入矩阵的每一列序列进行一维离散小波分解
     tmp2=x(:,k);
     [ca2,cd2]=mydwt(tmp2,lpd,hpd,1);
      x(:,k)=[ca2,cd2]; % 将分解所得系数存入矩阵x中，得到[LL,Hl;LH,HH]
end
LL=x(1:row/2,1:col/2); % LL是矩阵x的左上角部分
LH=x(row/2+1:row,1:col/2); % LH是矩阵x的左下角部分
HL=x(1:row/2,col/2+1:col); % HL是矩阵x的右上角部分
HH=x(row/2+1:row,col/2+1:col); % HH是矩阵x的右下角部分

function y = myidwt(cA,cD,lpr,hpr)
% 函数 MYIDWT() 对输入的小波分解系数进行逆离散小波变换，重构出信号序列 y
% 输入参数：cA —— 平均部分的小波分解系数；
% cD —— 细节部分的小波分解系数；
% lpr、hpr —— 重构所用的低通、高通滤波器。
 
lca=length(cA); % 求出平均、细节部分分解系数的长度
lcd=length(cD);
 
while (lcd)>=(lca) % 每一层重构中，cA 和 cD 的长度要相等，故每层重构后，
% 若lcd小于lca，则重构停止，这时的 cA 即为重构信号序列 y 。
upl=upspl(cA); % 对平均部分系数进行上抽样
cvl=conv(upl,lpr); % 低通卷积
 
cD_up=cD(lcd-lca+1:lcd); % 取出本层重构所需的细节部分系数，长度与本层平均部分系数的长度相等
uph=upspl(cD_up); % 对细节部分系数进行上抽样
cvh=conv(uph,hpr); % 高通卷积
 
cA=cvl+cvh; % 用本层重构的序列更新cA，以进行下一层重构
cD=cD(1:lcd-lca); % 舍弃本层重构用到的细节部分系数，更新cD
lca=length(cA); % 求出下一层重构所用的平均、细节部分系数的长度
lcd=length(cD);
end % lcd < lca，重构完成，结束循环
y=cA; % 输出的重构序列 y 等于重构完成后的平均部分系数序列 cA
 
function y=upspl(x)
% 函数 Y = UPSPL(X) 对输入的一维序列x进行上抽样，即对序列x每个元素之间
% 插零，例如 x=[x1,x2,x3,x4],上抽样后为 y=[x1,0,x2,0,x3,0,x4];
 
N=length(x); % 读取输入序列长度
M=2*N-1; % 输出序列的长度是输入序列长度的2倍再减一
for i=1:M % 输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列元素
if mod(i,2)
y(i)=x((i+1)/2);
else
y(i)=0;
end
end
 

function y=myidwt2(LL,HL,LH,HH)
% 基于Haar小波的二维小波分解和重构过程：
%函数 MYIDWT2() 对输入的子矩阵序列进行逆小波变换，重构出矩阵 y
% 输入参数：LL,HL,LH,HH —— 是四个大小均为 r*c 维的矩阵
% 输出参数：y —— 是一个大小为 2r*2c 维的矩阵
lpr=[1 1];hpr=[1 -1]; % 默认的低通、高通滤波器
tmp_mat=[LL,HL;LH,HH]; % 将输入的四个矩阵组合为一个矩阵
[row,col]=size(tmp_mat); % 求出组合矩阵的行列数
for k=1:col % 首先对组合矩阵tmp_mat的每一列，分开成上下两半
ca1=tmp_mat(1:row/2,k); % 分开的两部分分别作为平均系数序列ca1、细节系数序列cd1
cd1=tmp_mat(row/2+1:row,k);
tmp1=myidwt(ca1,cd1,lpr,hpr); % 重构序列
yt(:,k)=tmp1; % 将重构序列存入待输出矩阵 yt 的相应列，此时 y=[L|H]
end
 
for j=1:row % 将输出矩阵 y 的每一行，分开成左右两半
ca2=yt(j,1:col/2); % 分开的两部分分别作为平均系数序列ca2、细节系数序列cd2
cd2=yt(j,col/2+1:col);
tmp2=myidwt(ca2,cd2,lpr,hpr); % 重构序列
yt(j,:)=tmp2; % 将重构序列存入待输出矩阵 yt 的相应行，得到最终的输出矩阵 y=yt
end
y=yt;

function k=trancetion(img)
%将图像尺寸变为偶数
[m,n]=size(img);
if mod(m,2)==0
  w=m;
else
  w=m+1;    
end
if mod(n,2)==0
    w1=n;
else
    w1=n+1;
end
k=imresize(img,[w,w1]);

function Ioutput=MirrorImage(Iinput,wmax)
% to create a mirror image, wmax is the number of the extended rows or columns on every side of the image
% the Iinput should be double
[high,wid]=size(Iinput);
Ioutput=zeros(high+2*wmax,wid+2*wmax);
Ioutput(1+wmax:high+wmax,1+wmax:wid+wmax)=Iinput;
clear Iinput;
Ioutput(1:wmax,1+wmax:wid+wmax)=Ioutput(2*wmax:-1:1+wmax,1+wmax:wid+wmax);
Ioutput(high+wmax+1:high+2*wmax,1+wmax:wid+wmax)=Ioutput(high+wmax:-1:high+1,1+wmax:wid+wmax);
Ioutput(:,1:wmax)=Ioutput(:,2*wmax:-1:1+wmax);
Ioutput(:,wid+wmax+1:wid+2*wmax)=Ioutput(:,wid+wmax:-1:wid+1);