【转】小波图像分解与重构

原文地址：【转】小波图像分解与重构 作者：freedesert

小波图像融合综述（1）

分类： 小波变换 2007-11-14 11:15  4728人阅读  评论(6)  收藏  举报

       图像融合是将两幅或多幅图像融合在一起，以获取对同一场景的更为精确、更为全面、更为可靠的图像描述。融合算法应该充分利用各原图像的互补信息，使融合后的图像更适合人的视觉感受，适合进一步分析的需要；并且应该统一编码，压缩数据量，以便于传输。
图像融合可分为三个层次:
       1.    像素级融合
       2.    特征级融合
       3.    决策级融合
       其中像素级融合是最低层次的融合，也是后两级的基础。它是将各原图像中对应的像素进行融合处理，保留了尽可能多的图像信息， 精度比较高， 因而倍受人们的重视。像素级的图像融合方法大致可分为三大类:
       1.    简单的图像融合方法
       2.    基于塔形分解（如Laplace塔形分解、比率塔等）的图像融合方法
       3.    基于小波变换的图像融合方法
       小波变换是图像的多尺度、多分辨率分解，它可以聚焦到图像的任意细节，被称为数学上的显微镜。近年来，随着小波理论及其应用的发展，已将小波多分辨率分解用于像素级图像融合。小波变换的固有特性使其在图像处理中有如下优点：
       1.    完善的重构能力，保证信号在分解过程中没有信息损失和冗余信息；
       2.    把图像分解成平均图像和细节图像的组合，分别代表了图像的不同结构，因此容易提取原始图像的结构信息和细节信息；
       3.    具有快速算法，它在小波变换中的作用相当于FFT算法在傅立叶变换中的作用，为小波变换应用提供了必要的手段；
       4.    二维小波分析提供了与人类视觉系统方向相吻合的选择性图像。

——像素级图像融合的主要步骤
       以两幅图像的融合为例。设A，B为两幅原始图像，F为融合后的图像。若对二维图像进行N层的小波分解，最终将有(3N+1)个不同频带，其中包含3N 个高频子图像和1个低频子图像。其融合处理的基本步骤如下： 
       （1）对每一原图像分别进行小波变换，建立图像的小波塔型分解； 
       （2）对各分解层分别进行融合处理。各分解层上的不同频率分量可采用不同的融合算子进行融合处理，最终得到融合后的小波金字塔； 
       （3）对融合后所得小波金字塔进行小波重构，所得到的重构图像即为融合图像。

                                                                                   图 1 

      在图像融合过程中，小波基的种类和小波分解的层数对融合效果有很大的影响，对特定的图像来说，哪一种小波基的融合效果最好，分解到哪一层最合适，都是需要考虑的问题。为此可以通过引入融合效果的评价来构成一个闭环系统。如图2所示。

自己动手编写小波信号分解与重构的Matlab程序

分类： 小波变换 2007-11-13 10:56  9744人阅读  评论(31)  收藏  举报

 

P.S.：（2008-09-01）感谢网友‘李明杨艳’指出了本文程序中一维信号小波分解重构程序mydwt和myidwt存在的一个大Bug，现已修正，请参见今天发表的文章《一维信号的小波分解与重构程序》。

P.S.：（2008-06-05）去年11月发布了一系列有关小波变换和图像处理的文章，把学习小波过程中的心得体会和编写的程序放在网上和大家共享交流。半年来，感谢大家的关注和帮助，在相互的讨论交流中，我不断地从大家提出的问题中拓展自己的知识面，对小波的理论及其应用有了更深入的了解和掌握。根据和大家讨论交流中发现的问题，对博客中的程序进行修正。有关小波图像分解和重构的两篇文章中分享的程序，存在下列问题：

（1）程序所用的小波函数只有非标准的Haar小波，其滤波器组为 Lo_D=[1/2 1/2], Hi_D=[-1/2 1/2]，是固化在 mydwt2.m 的程序中的，不能选择其他的小波函数；

（2）非标准的Haar小波，其分解出来的系数矩阵中，高频系数的细节内容（轮廓、边缘等特征）不明显；

（3）函数 mydwt2 中列变换的矩阵对象为输入矩阵，这是错误的，其矩阵对象应该是行变换后的缓存矩阵；

（4）函数 mydwt2 的输出用[LL,HL,LH,HH]表示，不是很规范，应改为[cA,cV,cH,cD]来表示，即一级小波变换输出的系数矩阵有4个部分：平均部分、垂直细节部分、水平细节部分和对角线细节部分。

（5）函数 mywavedec2 的输出 y 是与输入矩阵 x 相同大小的矩阵，并且已将N级分解后所有的平均、细节系数组合成一体的。实际上，这种定义只对Haar小波有效。

（6）原程序中要调用 modmat 函数对图像矩阵进行修剪，使之能被 2 的 N 次方整除，主要是为了生成塔式结构图像而设的，对上述问题修正后，这个 modmat 函数已不需使用了。

针对上述问题，我对程序作了修正，发布在今天的3篇文章里，请大家点击查看。新修正的程序更为简洁易懂，功能也有所增强，可以用任意的小波函数进行小波分解，可根据小波分解系数矩阵重构出指定分解级的低频系数和原始图像。

         1、《小波图像分解与重构程序存在的问题与解决办法》
http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2513126.aspx

2、《小波图像分解 Matlab 程序 - V2.0版》
http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2513865.aspx

3、《小波图像重构 Matlab 程序 - V2.0版》
http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2514119.aspx

 

（2007-11-13）Matlab小波分析工具箱丰富的函数和强大的仿真功能为我们学习小波、用好小波提供了方便、快捷的途径，但是，如果我们要深入掌握小波分析的原理，真正学好、用好小波，就应该尽量用自己编写的程序去实现小波变换和信号分析，尽量在自己的程序中少调用Matlab提供的函数，多用自己的理解去编写相关的小波函数，这样的过程是一个探索、求知的过程，更能让我们体会到小波的强大和学习的乐趣。下面，我把自己编写的小波一维、二维信号分解和重构Matlab程序共享出来，也希望有朋友共享自编的程序，共同学习，提高程序的效率和简洁性。

 

首先要说明的一点是，虽然是自己编写Matlab程序，但并不是说一点也不用Matlab的自带函数。我们要编写的是实现小波变换的主要功能函数，而绘图等基本功能还是要用到Matlab函数的。而且，根据小波变换的滤波器组原理，原始信号要通过低通、高通滤波器处理，这里就涉及到卷积这一运算步骤。卷积——FFT算法的实现，相信很多朋友都能用Matlab、C语言等来实现，不过与Matlab自带的用机器语言编写的FFT程序相比，运算速度一般会慢几倍、几十倍。所以，我的程序里边涉及卷积的就直接调用Matlab的conv()函数了。

 

我们知道，小波变换的一级分解过程是，原始信号分别进行低通、高通滤波，再分别进行二元下抽样，就得到低频、高频（也称为平均、细节）两部分系数；而多级分解则是对上一级分解得到的低频系数再进行小波分解，是一个递归过程。以下是一维小波分解的程序：

 

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 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

   

    我们知道，二维小波分解重构可以用一系列的一维小波分解重构来实现。以下程序是基于Haar小波的二维小波分解和重构过程：

 

function [LL,HL,LH,HH]=mydwt2(x);
% 函数 MYDWT2() 对输入的r*c维矩阵 x 进行二维小波分解，输出四个分解系数子矩阵[LL,HL,LH,HH]
% 输入参数：x —— 输入矩阵，为r*c维矩阵。
% 输出参数：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=myidwt2(LL,HL,LH,HH);
% 函数 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;

   

    我编写程序时尽量实现模块化，方便不同功能程序之间的调用。以上程序实现了一维、二维小波分解与重构的功能，接下来我们就可以此为基础，探讨小波信号压缩和图像融合等技术了。

用自编的程序实现小波图像分解与重构

分类： 小波变换 2007-11-13 12:08  7634人阅读  评论(25)  收藏  举报

去年11月发布了一系列有关小波变换和图像处理的文章，把学习小波过程中的心得体会和编写的程序放在网上和大家共享交流。半年来，感谢大家的关注和帮助，在相互的讨论交流中，我不断地从大家提出的问题中拓展自己的知识面，对小波的理论及其应用有了更深入的了解和掌握。根据和大家讨论交流中发现的问题，对博客中的程序进行修正。有关小波图像分解和重构的两篇文章中分享的程序，存在下列问题：

（1）程序所用的小波函数只有非标准的Haar小波，其滤波器组为 Lo_D=[1/2 1/2], Hi_D=[-1/2 1/2]，是固化在 mydwt2.m 的程序中的，不能选择其他的小波函数；

（2）非标准的Haar小波，其分解出来的系数矩阵中，高频系数的细节内容（轮廓、边缘等特征）不明显；

（3）函数 mydwt2 中列变换的矩阵对象为输入矩阵，这是错误的，其矩阵对象应该是行变换后的缓存矩阵；

（4）函数 mydwt2 的输出用[LL,HL,LH,HH]表示，不是很规范，应改为[cA,cV,cH,cD]来表示，即一级小波变换输出的系数矩阵有4个部分：平均部分、垂直细节部分、水平细节部分和对角线细节部分。

（5）函数 mywavedec2 的输出 y 是与输入矩阵 x 相同大小的矩阵，并且已将N级分解后所有的平均、细节系数组合成一体的。实际上，这种定义只对Haar小波有效。

（6）原程序中要调用 modmat 函数对图像矩阵进行修剪，使之能被 2 的 N 次方整除，主要是为了生成塔式结构图像而设的，对上述问题修正后，这个 modmat 函数已不需使用了。

针对上述问题，我对程序作了修正，发布在今天的3篇文章里，请大家点击查看。新修正的程序更为简洁易懂，功能也有所增强，可以用任意的小波函数进行小波分解，可根据小波分解系数矩阵重构出指定分解级的低频系数和原始图像。

1、《小波图像分解与重构程序存在的问题与解决办法》
http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2513126.aspx

2、《小波图像分解 Matlab 程序 - V2.0版》
http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2513865.aspx

3、《小波图像重构 Matlab 程序 - V2.0版》
http://blog.csdn.net/chenyusiyuan/archive/2008/06/05/2514119.aspx

 

      上一篇文章中我们实现了小波的一维、二维信号分解与重构，其中的二维信号分解与重构，只要稍作修改，就可以实现图像的分解和重构了。修改的工作，主要是对图像信号进行规范化处理、数据格式转换和绘图细节处理等。

 

       简单起见，我们从黑白（灰度）图像的分解、重构说起，因为彩色图像的处理要复杂一点。在本文中，我们使用著名的Lena图作为原始图像。

图1

首先，为了实现图像的N层分解，对一幅m行n列的黑白图像，我们要对其进行规范化处理，使其能被2的N次方整除。以下的modmat() 函数实现此功能：

 

function y=modmat(x,dim)
% 函数 MODMAT() 对输入矩阵x进行规范化，使其行列数均能被 2^dim 整除
% 输入参数：x —— r*c 维矩阵；
%           dim —— 矩阵重构的维数
% 输出参数：y —— rt*ct 维矩阵，mod(rt,2^dim)=0，mod(ct,2^dim)=0

[row,col]=size(x);          % 求出输入矩阵的行列数row,col
rt=row - mod(row,2^dim);    % 将row,col分别减去本身模 2^dim 得到的数
ct=col - mod(col,2^dim);    % 所得的差为rt、ct，均能被 2^dim 整除
y=x(1:rt,1:ct);             % 输出矩阵 y 为输入矩阵 x 的 rt*ct 维子矩阵

 

然后，将规范化后的图像的数据格式由适合显示图像的uint8格式转换为适合数值处理的double格式，再调用二维小波分解函数进行图像分解，最后为了清晰地显示分解图像的塔式结构，在图像的相应区域绘制若干分界线。具体程序如下：

 

function y=mywavedec2(x,dim)
% 函数 MYWAVEDEC2() 对输入矩阵 x 进行 dim 层分解，得到相应的分解系数矩阵 y
% 输入参数：x —— 输入矩阵；
%           dim —— 分解层数。
% 输出参数：y —— 分解系数矩阵。

x=modmat(x,dim);            % 首先规范化输入矩阵，使其行列数均能被 2^dim 整除

subplot(121);imshow(x);title('原始图像');   % 画出规范化后的源图像
[m,n]=size(x);              % 求出规范化矩阵x的行列数
xd=double(x);               % 将矩阵x的数据格式转换为适合数值处理的double格式

for i=1:dim
    xd=modmat(xd,1);
    [dLL,dHL,dLH,dHH]=mydwt2(xd);   % 矩阵小波分解
    tmp=[dLL,dHL;dLH,dHH];          % 将分解系数存入缓存矩阵
    xd=dLL;                         % 将缓存矩阵左上角部分的子矩阵作为下一层分解的源矩阵
    [row,col]=size(tmp);            % 求出缓存矩阵的行列数
    y(1:row,1:col)=tmp;             % 将缓存矩阵存入输出矩阵的相应行列
end

yd=uint8(y);            % 将输出矩阵的数据格式转换为适合显示图像的uint8格式
for i=1:dim             % 对矩阵 yd 进行分界线处理，画出分解图像的分界线
    m=m-mod(m,2);
    n=n-mod(n,2);
    yd(m/2,1:n)=255;
    yd(1:m,n/2)=255;
    m=m/2;n=n/2;
end
subplot(122);imshow(yd);title([ num2str(dim) ' 维小波分解图像']);

 

以下是程序的运行结果。

图2

       上述的图像分解程序，其输出数据是double格式的，以便作为重构程序的输入。

接下来我们讨论图像的重构。我编写的重构程序中，为了比较分解图像和重构图像，首先绘出经过小波分解的图像，然后再进行重构。在绘制分解图像和重构图像的过程中，要注意数据格式的转换。

 
function y=mywaverec2(x,dim)

% 函数 MYWAVEREC2() 对输入的分解系数矩阵x进行 dim 层重构，得到重构矩阵 y

% 输入参数：x ——分解系数矩阵；

%           dim ——重构层数；

% 输出参数：y ——重构矩阵。

% 绘制分解图像

xd=uint8(x);            % 将输入矩阵的数据格式转换为适合显示图像的uint8格式

[m,n]=size(x);          % 求出输入矩阵的行列数

for i=1:dim             % 对转换矩阵xd进行分界线处理

    m=m-mod(m,2);

    n=n-mod(n,2);

    xd(m/2,1:n)=255;

    xd(1:m,n/2)=255;

    m=m/2;n=n/2;

end

figure;

subplot(121);imshow(xd);title([ num2str(dim) ' 层小波分解图像']); % 画出带有分界线的分解图像

% 重构图像

xr=double(x);           % 将输入矩阵的数据格式转换回适合数值处理的double格式

[row,col]=size(xr);     % 求出转换矩阵xr的行列数

for i=dim:-1:1          % 重构次序是从内层往外层进行，所以先抽取矩阵 xr 的最内层分解矩阵进行重构

    tmp=xr(1:floor(row/2^(i-1)),1:floor(col/2^(i-1)));       % 重构的内层矩阵的行列数均为矩阵xr的2^(i-1)

    [rt1,ct1]=size(tmp);                         % 读取待重构矩阵 tmp 的行列数

    rt=rt1-mod(rt1,2);ct=ct1-mod(ct1,2);

    rLL=tmp(1:rt/2,1:ct/2);                    % 将待重构矩阵 tmp 分解为四个部分

    rHL=tmp(1:rt/2,ct/2+1:ct);

    rLH=tmp(rt/2+1:rt,1:ct/2);

    rHH=tmp(rt/2+1:rt,ct/2+1:ct);

    tmp(1:rt,1:ct)=myidwt2(rLL,rHL,rLH,rHH);              % 将重构结果返回到矩阵 tmp

    xr(1:rt1,1:ct1)=tmp;       % 把矩阵 tmp 的数据返回到矩阵 xr 的相应区域，准备下一个外层的重构

end

y=xr;                    % 重构结束后得到的矩阵xr即为输出矩阵 y

yu=uint8(xr);            % 将矩阵xr的数据格式转换为适合显示图像的uint8格式

subplot(122);imshow(yu);title('小波重构图像 ');

图3