多尺度小波分解得到的系数进行处理之后如何进行重构的问题说明

来源:互联网 发布:linux怎么安装make 编辑:程序博客网 时间:2024/04/29 03:17

毕设时候做的题目是《基于分数傅里叶变换和小波变换的图像加密解密》,当时分数傅里叶变换只是简单的应用了下,没留下什么太大的印象。反而是在应用小波的一些处理过程中,因为先后出现了很多问题,留了一些印象。记的最清楚的是一个这样的问题:试图使用二维多尺度离散小波变换,然后对得到的小波系数进行Arnold置乱处理,如果利用matlab子带的wavedec2,可以很容易实现分解(格式为:[c,s]=wavedec2(X,N,'wname')或者[c,s]=wavedec2(X,N,Lo_D,Hi_D)),然后利用appcoef2和detcoef2可以很容易地将低频和高频提取出来,例如:提取低频系数可以有A = appcoef2(C,S,'wname',N),提取水平方向的高频系数可以有H = detcoef2('h',C,S,N),提取垂直方向高频系数可以有
V = detcoef2('v',C,S,N);,提取对角线方向的高频系数可以有D= detcoef2('d',C,S,N),其中N为分解尺度。但是对提取出来的系数进行处理后,然后如何畸形重构呢?Matlab提供waverec2函数进行重构,格式为:X=waverec2(C,S,'wname') 或者X=waverec2(C,S,Lo_R,Hi_R) 。但是,处理后的系数是如何构成C和S呢?

当然可以使用维单尺度离散小波变换,通过迭代进行对低频进行单尺度离散小波分解,也可以达到对原图像进行多尺度分解的目的。而且若使用二维单尺度离散小波变换,各个系数的提取以及单支重构都有matlab现成的函数,很容易进行操作。唯一存在的问题是当避开二维多尺度离散变换函数,试图利用二维单尺度离散小波变换函数dwt2对每次分解后得到的低频系数不断进行处理时,在某些小波基的情况下会造成重构的失败。当时,也不知道有imresize这个函数,选的图像集就全部选了256*256大小的图片。汗。

还在各种QQ群里骚扰各路大神,有个同学给出了matlab自带的help里面的截图,说了句c里面存的就是你要的东西啊。额,不给自己找借口了,反正结果就是,最终上面的问题还是被自己堂而皇之的规避了过去。现在又遇上了,就猛的想起来哪位同学的话,就想试试。也正好好好地领略一下matlab的help的强大。当然,还是借助了网上的一些解释(一并谢过了)总算是弄得明白了很多。不怕被笑话,也总结一下。

在Mtlab里面help下wavedec2函数出来的说明里面,最重要的是下面的部分:

 The output wavelet 2-D decomposition structure [C,S] contains the wavelet decomposition vector C and the corresponding bookeeping matrix S.  Vector C is organized as:
C = [  A(N) |  H(N)  |  V(N)  |  D(N) |  H(N-1) | V(N-1) | D(N-1) | ...  | H(1) | V(1) | D(1) ]. where A, H, V, D, are row vectors such that:  A = approximation coefficients,  H = hori. detail coefficients, V = vert. detail coefficients, D = diag. detail coefficients, each vector is the vector column-wise storage of a matrix. Matrix S is such that:  S(1,:) = size of app. coef.(N) , S(i,:) = size of det. coef.(N-i+2) for i = 2,...,N+1 and S(N+2,:) = size(X).

这段话说的是,输出的二维小波分解结构为[C, S],其中C为各层分解系数,S为各层分解系数的长度,即大小。其中,C的结构为[  A(N) |  H(N)  |  V(N)  |  D(N) |  H(N-1) | V(N-1) | D(N-1) | ...  | H(1) | V(1) | D(1) ],A(N)代表第N层低频系数,H(N),V(N)和D(N)分别代表第N层的水平高频系数、垂直高频系数和对角高频系数,类推到H(1),V(1)和D(1)。C是一个行向量,由这些向量做成,而每个向量是一个矩阵的每列转置的组合存储(column-wise storage of a matrix)。如果我们知道每个向量的大小,然后恢复成矩阵,就可以得到相应的分解系数。向量的大小可以有S得到。

S储存的是各层分解系数的长度,确切的说是原矩阵的行数和列数。具体来说,第一行是A(N)的行数和列数,第二行是H(N)的大小,因为同一层尺度的分解系数的大小相同,所以,第二行也是V(N)或者D(N)的大小,依次类推,直到H(1)(或者V(1),D(1))的大小,注意,还有最后一行为X的长度,即输入图像的大小。于是,S中存储的行数的大小和列数的大小相乘得到的就是该层分解得到的储存在C中的向量的大小。这样,根据S可以求出C中相应的向量的大小。

由上面可以想到,对提取出来的各层系数(为矩阵形式)进行处理之后,利用reshape函数或者直接转换为行向量,然后按照同样的顺序就可以形成处理后的C,S不变,然后利用matlab的waverec2函数就可以进行图像的重构了。

下面,我把自己对文章《Change Detection in Synthetic Aperture Radar Images based on Image Fusion and Fuzzy Clustering》的仿真中小波融合的代码部分贴出来,算是给出一个实例吧。

(友情提示:可直接看最后重构部分的示例)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 本例是针对两幅大小相同的灰度图像,分别对低频系数和高频系数进行融合

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;
img1=imread('img1.bmp');img2=imread('img2.bmp'); %读入需要融合的两幅图像
N=3;wname='bior3.7'; %小波分解级数及类型
[C1,S1]=wavedec2(img1,3,wname);
[C2,S2]=wavedec2(img2,3,wname); %对两幅图像进行多尺度离散小波分解
%分别对每幅图像提取各尺度的各个系数
%这里指,第三尺度的低频系数,以及各个尺度的高频系数,包括水平、垂直和对角线方向
%尺度1:提取高频系数
cH11=detcoef2('h',C1,S1,1);cV11=detcoef2('v',C1,S1,1);cD11=detcoef2('d',C1,S1,1);
cH12=detcoef2('h',C2,S2,1);cV12=detcoef2('v',C2,S2,1);cD12=detcoef2('d',C2,S2,1);
%对系数进行处理,规则见原文章
[m1,n1]=size(cH11); %当然,也可以通过S得到其大小,例如m1=S1(4,1),n1=S1(4,2);
%用3*3滑动窗口求每个象素的能量
H11=nlfilter(cH11,[z1 z1],@myfunction3);
V11=nlfilter(cV11,[z1 z1],@myfunction3);
D11=nlfilter(cD11,[z1 z1],@myfunction3);
H12=nlfilter(cH12,[z1 z1],@myfunction3);
V12=nlfilter(cV12,[z1 z1],@myfunction3);
D12=nlfilter(cD12,[z1 z1],@myfunction3);
for i=1:m1
    for j=1:n1         
          if abs(H11(i,j))<abs(H12(i,j))
                cH12(i,j)=cH11(i,j);
          else
                cH12(i,j)=cH12(i,j);
          end
          if abs(V11(i,j))<abs(V12(i,j))
                cV12(i,j)=cV11(i,j);
          else
                cV12(i,j)=cV12(i,j);
          end
          if abs(D11(i,j))<abs(D12(i,j))
                cD12(i,j)=cD11(i,j);
          else
                cD12(i,j)=cD12(i,j);
          end
    end
end
%尺度2:提取高频系数
cH21=detcoef2('h',C1,S1,2);cV21=detcoef2('v',C1,S1,2);cD21=detcoef2('d',C1,S1,2);
cH22=detcoef2('h',C2,S2,2);cV22=detcoef2('v',C2,S2,2);cD22=detcoef2('d',C2,S2,2);
%对系数进行处理
[m2,n2]=size(cH21); %同样,也可以通过S得到其大小
%用3*3滑动窗口求每个象素的能量
H21=nlfilter(cH21,[z1 z1],@myfunction3);
V21=nlfilter(cV21,[z1 z1],@myfunction3);
D21=nlfilter(cD21,[z1 z1],@myfunction3);
H22=nlfilter(cH22,[z1 z1],@myfunction3);
V22=nlfilter(cV22,[z1 z1],@myfunction3);
D22=nlfilter(cD22,[z1 z1],@myfunction3);
for i=1:m2
    for j=1:n2
          if abs(H21(i,j))<abs(H22(i,j))
              cH22(i,j)=cH21(i,j);
          else
              cH22(i,j)=cH22(i,j);
          end
          if abs(V21(i,j))<abs(V22(i,j))
              cV22(i,j)=cV21(i,j);
          else
              cV22(i,j)=cV22(i,j);
          end
          if abs(D21(i,j))<abs(D22(i,j))
              cD22(i,j)=cD21(i,j);
          else
              cD22(i,j)=cD22(i,j);
          end
    end
end
%尺度3:提取低频与高频系数
cA31=appcoef2(C1,S1,wname,3); 提取尺度3的低频系数
cH31=detcoef2('h',C1,S1,3);cV31=detcoef2('v',C1,S1,3);cD31=detcoef2('d',C1,S1,3);
cA32=appcoef2(C2,S2,wname,3);提取尺度3的低频系数
cH32=detcoef2('h',C2,S2,3);cV32=detcoef2('v',C2,S2,3);cD32=detcoef2('d',C2,S2,3);
%对系数进行处理
[m3,n3]=size(cH31); %当然,也可以通过S得到其大小
%用3*3滑动窗口求每个象素的能量
H31=nlfilter(cH31,[z1 z1],@myfunction3);
V31=nlfilter(cV31,[z1 z1],@myfunction3);
D31=nlfilter(cD31,[z1 z1],@myfunction3);
H32=nlfilter(cH32,[z1 z1],@myfunction3);
V32=nlfilter(cV32,[z1 z1],@myfunction3);
D32=nlfilter(cD32,[z1 z1],@myfunction3);
A31=nlfilter(cA31,[z1 z1],@myfunction3);
A32=nlfilter(cA32,[z1 z1],@myfunction3);
for i=1:m3
    for j=1:n3
          if abs(H31(i,j))<abs(H32(i,j))
               cH32(i,j)=cH31(i,j);
          else
               cH32(i,j)=cH32(i,j);
          end
          if abs(V31(i,j))<abs(V32(i,j))
               cV32(i,j)=cV31(i,j);
          else
               cV32(i,j)=cV32(i,j);
          end
          if abs(D31(i,j))<abs(D32(i,j))
               cD32(i,j)=cD31(i,j);
          else
               cD32(i,j)=cD32(i,j);
          end
          if abs(A31(i,j))<abs(A32(i,j))
               cA32(i,j)=(cA31(i,j)+cA32(i,j))/2;
          else
               cA32(i,j)=(cA31(i,j)+cA32(i,j))/2;                       
          end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% ------------------------------------------------------ 重构--------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

len1=S1(1,1)*S1(1,2); % A(3),H(3),V(3),D(3)的大小,当然也可以用之前得到的m1,n1等
len2=S1(3,1)*S1(3,2); % H(2),V(2),D(2)的大小
len3=S1(4,1)*S1(4,2); % H(1),V(1),D(1)的大小
C(1:len1)=cA32(1:end); % 将A(3)放入相应的位置
C(len1+1:2*len1)=cH32(1:end); % 将H(3),V(3),D(3)放入相应的位置
C(2*len1+1:3*len1)=cV32(1:end);
C(3*len1+1:4*len1)=cD32(1:end);
C(4*len1+1:4*len1+len2)=cH22(1:end); % 将H(2),V(2),D(2)放入相应的位置
C(4*len1+len2+1:4*len1+2*len2)=cV22(1:end);
C(4*len1+2*len2+1:4*len1+3*len2)=cD22(1:end);
C(4*len1+3*len2+1:4*len1+3*len2+len3)=cH12(1:end); % 将H(1),V(1),D(1)放入相应的位置
C(4*len1+3*len2+len3+1:4*len1+3*len2+2*len3)=cV12(1:end);

C(4*len1+3*len2+2*len3+1:4*len1+3*len2+3*len3)=cD12(1:end);

% 得到系数处理之后构成的新的C

IMG=waverec2(C,S1,wname);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% ------------------------------------------------------ END--------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0 0
原创粉丝点击