使用MATLAB设计小波变换程序中的若干问题

suffer · 发表于 2005-8-11 20:35

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

<TABLE height="85%" width="95%" align=center border=0>

<TR>
<TD vAlign=top width=* height="100%">在使用MATLAB完成小波变换程序和通过阈值来压缩图像的过程中，我和许多同学都是边学边用，是从一个接一个的问题中逐步理解小波和MATLAB编写程序的。因此我愿意就个人遇到和解决问题的经验与大家讨论，希望能够对遇到同样问题的人有所帮助。在清华大学林福宗老师倡导的网上互动的学习方式中，老师同学的开诚布公的讨论，尤其是林老师启发大家对出现‘问题’采取的态度和做法，对我今后成长为一名合格的清华的研究生意义重大，谨以此文表示对他指导关心的敬意！ 内容简介 本文分为三部分: 如何使用MATLAB设计小波标准与标准分解；如何完成使用小波变换压缩图像；仍需探讨的问题。每一部分主要以问题和例子的形式讲述，为了便于参照，附录部分给出部分源代码供大家参考指正。 我个人在完成作业的时候，走了许多弯路，最后，反复读老师的第3章讲义的第3--5节，明白了小波变换的一些非常基础性的知识，正是如此，我主要是按照第三章的例子，做出使用MATLAB进行Haar变换的实例，这至少对完成非标准Haar小波的3级分解与合成没有任何问题，而且计算速度远远比使用一维变换快，甚至超过dwt2和wavedec2。再者，是通过这些简短的例子的学习，也能从实践的角度理解小波变换的基础知识和道理，我觉得掌握小波的知识要比简单使用dwt2直接分析出结果更重要。为了让使用一维变换或打算采用卷积完成分解重构程序的同学有所参考，也给出了我以前写的代码和设想供有兴趣者参考。在图像压缩的任务中，我个人建议采用wdencmp，原因是简单可*，能够生成老师要求小波压缩与重构的演示及PNG文件。在最后部分着重对使用小波压缩之后重构图像的PNG文件会变大问题进行了简单的猜测性的解释和讨论，希望感兴趣的同学能深入研究。 如何使用MATLAB设计小波标准与标准分解 使用Haar和Db9小波编写图像文件的标准和非标准分解重构程序 实质上，任务书中规定首先要编写４种程序之一：哈尔（haar）小波的标准分解重构程序；Daubechies 9小波的标准分解重构程序；哈尔（haar）小波的非标准分解重构程序；Daubechies 9小波的非标准分解重构程序。如果采用标准方法，需要生成图３－２５的系列图像。如果采用非标准方法则需要生成图３－２６的系列图像。 我个人编写了非标准的分解重构程序，并总结网上同学们公开的方法，认为可以有３种方法实现。但前提是必须仔细阅读林老师推荐的补充教材第三章“小波与小波变换”的3.3、3.4和3.5节（１４－２７页）。如果能使用MATLAB简单地实践一下教材的例子，至少完成哈尔小波的标准与非标准程序相当简单。为此我主要介绍如何实现老师教材给我们的例子是如何在MATLAB中实现的。然后，再引申到使用MATLAB的函数实现方法。为了描述问题简便，使用黑体表示需在MATLAB的命令窗口（command window）的输入部分。如果不熟悉下面出现的函数功能和使用方法，请再命令窗口中使用 help 函数名，如help dwt2，能够得到英文的功能和用法说明。 哈尔小波变换的MATLAB实例 先介绍MATLAB的矩阵(Matrix)和向量(Vector)的赋值方法：在command window的>>符号下输入： 一维向量：　Ｉ = [ 9 7 3 5] 　　　　　　Ｆ = [2, 5, 8, 9, 7, 4, -1, 1] 二维矩阵： A = [ 64 2 3 61 60 6 7 57 9 55 54 12 13 51 50 16 17 47 46 20 21 43 42 24 40 26 27 37 36 30 31 33 32 34 35 29 28 38 39 25 41 23 22 44 45 19 18 48 49 15 14 52 53 11 10 56 8 58 59 5 4 62 63 1 ] 　 Ｉ和Ｆ现在就是［例3.1］［例3.2］一维向量,　A就是3.5.1中的图像矩阵。在进行Haar小波的一维与二维变换（分解）之前，首先介绍一下MATLAB的矩阵运算。 MATLAB是Matrix Laboratory 的合写，它对矩阵运算之功能堪称一流。由于使用矩阵描述问题更象数学表达式，所以编写的程序不仅高效，更易读。所以MATLAB程序应该尽量使用矩阵直接描述。如 M是一4X4的矩阵，则B=I*M就完成了使用两个for循环编写乘法的程序。 分析老师的［例3.1］，那么可以得到如下结果： [（9+7）*1/2 （3+5）* 1/2 （9-7）*1/2 （3-5）* 1/2] 如果把其看作矩阵方式的乘法，那么令 M为其表述求取平均和差值的系数矩阵，则输入： M = [ 1/2 0 1/2 0 1/2 0 -1/2 0 0 1/2 0 1/2 0 1/2 0 -1/2 ] 把上述的M输入到MATLAB中,然后使用： C=I*M 得到的结果就是 [ 8 4 1 -1]。这就是非规范化Haar小波的第一级分解系数。注意C的前两项[ 8 4 ] 就是近似系数（Approximation Coefficients）, 后两项[1 _1]就是细节系数（Detail Coefficients）. M就是Haar小波的非规范化（non-normalization）系数矩阵。如果我们执行： M*M_ (M_代表的是M的转置矩阵)，你会发现得到对角线为0.5其它为0的4X4矩阵，所以如果让它变成单位矩阵（对角线为1），必须把原来的1/2增大，变为1/sqrt(2), 所以执行： N=M*sqrt(2) 后得到新的系数矩阵 [ 0.7071 0 0.7071 0 0.7071 0 -0.7071 0 0 0.7071 0 0.7071 0 0.7071 0 -0.7071 ] 再执行 N * N_ 就会得到单位矩阵。N就是Haar小波的规范化（normalization）系数矩阵。现在执行: C*M_ 将得到 [4.5 3.5 1.5 2.5]，显然是因为M*M_的对角矩阵的值为0.5。现在使用 Cn=I*N 得到 [11.3137 5.6569 1.4142 -1.4142 ] 规范化的一维Haar小波一级分解系数。 I1=Cn*N_ 得到 [9 3 7 7] 规范化的一维Haar小波逆变换（重构）结果。可见得到这样规范化矩阵的好处是为了逆变换时不需要额外的运算，只是需要分解系数矩阵乘以规范化系数矩阵的转置。</TD></TR>
<TR>
<TD></TD></TR></TABLE>

suffer · 发表于 2005-8-11 20:35

如果需要对分解系数矩阵进行第2级分解，那么此时只对C或Cn的前2项的低频系数进行，此时的Haar的非规范化与规范化系数矩阵M1和N1分别为： M1=[ 1/2 1/2 1/2 &#0;1/2 ] N1=[ 1/2 1/2 1/2 &#0;1/2 ] 那么第2级分解可以用以下方式进行： C1=C(1:2)*M1 Cn1=Cn(1:2)*N1 结果是C1=[6　2 ],　Cn1= [12.0 4.0]。注意, C1(1:2)表示取出C中第１到第２个向量，符号“：”的含义是从１到２的意思，这是MATLAB截取数据的方法。逆变换： [C1*M1’*2　C(3:4)]*M’*2 [Cn1*N1’　Cn(3:4)]*N’ 这里采用了中括号［］的赋值方式直接将第２级的重构合并到第１级的系数中再重构原始向量I。 从上面的实验可以看出，除了输入哈尔小波系数矩阵M和N比较麻烦，整个变换和重构过程异常简单,这就启示我们只要编写一个Haar系数生成矩阵的函数，整个的分解问题即将化简。 问题1：如何生成Haar小波系数矩阵？ MATLAB的函数可以如以下形势构造：在MATLAB中使用File? New? M-file,输入 function [M]= haarmatrix (rows,cols, flag) %矩阵的行数 rows, 列数cols, 非规范化flag=０否则就是规范化矩阵 % M为函数返回的Haar系数矩阵 M=zeros(rows,cols);　　　　　　%首先制作rows X cols 的０矩阵 if flag== 0 　　　　　　　　 sv=0.5; %non-normalized else sv=1.0/sqrt(2.0); % normalized end half=rows/(2);　％矩阵从中间分开，左半部为近似系数部分，右半部为细节系数部分 for i= 1:2:half*2 M(i,ceil(i/2)) =sv; M(i+1,ceil(i/2))=sv; M(i,ceil(i/2)+half) =sv; M(i+1,ceil(i/2)+half)=-sv; end M=sparse(M) %生成稀疏矩阵，提高运算速度 注意分号的作用是不在command窗口上显示M的值。在输入后，请保存与函数名同名的文件名haarmatrix.m。现在来做［例3.2］： Ｆ = [2, 5, 8, 9, 7, 4, -1, 1] N= haarmatrix(8,8,1) C1=F*N N1= haarmatrix(4,4,1) C2= C1(1:4)*N1 N2= haarmatrix(2,2,1) C3= C2(1:2)*N2 按图3-21小波分解的层次结构合成最终系数 C=[C3, C2(3:4), C1(5:8)] 得到的C就是老师让我们使用第三章17页伪代码编写的一维小波变换分解的结果。其逆变换（重构）代码如下： Fr=[ [C3*N2' ,C2(3:4)]*N1' ,C1(5:8) ]* N' 注意理解最内部的中扩号实质上是使用C3逆变换的结果（低频）与C2的后２位细节系数一起逆变换得倒C1的前４位的低频近似系数，然后与C1的后４位细节系数构造出原始向量Fr。 为了验证我们的一维Haar变换程序，我们使用MATLAB提供的一维离散小波变换函数dwt来进行同样的计算： 第一级分解 [ca1, cd1]=dwt(F, ' haar' ) 第２级分解 [ca2,cd2]=dwt(ca1, ' haar' ) 第３级分解 [ca3,cd3]=dwt(ca2, ' haar' ) 按图3-21小波分解的层次结构合成最终系数 c=[ca3 cd3 cd2 cd1] 或直接使用MATLAB提供的１维多级分解程序 c=wavedec(F,3,'haar') 比较c与C　(注意MATLAB中大小写变量的区别)，可知我们成功构造了Haar小波分解合成程序。也许这是你可能觉得刚才的系数矩阵比起MATLAB提供的函数来说还是复杂了，但如果把上述过程函数化，做一个matrixDWT和matrixIDWT函数，使用一样简单，而且关键在于，我们使用矩阵方式解决２维图像的小波非标准变换的演示过程来说会带来更大便利。 问题2 如何进行图像的1级非标准分解和重构？ 如果我们利用2.2节开始时已经输入的18页的图像矩阵A，那么它的第一级二维变换的非标准分解过程只需３步：(为了方便核对MATLAB的dwt2分解结果，使用规范化系数矩阵) M1=haarmatrix(8,8,1) ; AR1=A*M1; 　 %这就是进逐行分解的图像 AC1=M1'*AR1; 　 %这就是进逐列分解的图像，也是一次２维变换的系数矩阵 使用MATLAB的dwt2来验证 [ca ch cv cd]=dwt2(A,'haar') c=[ca cv ; ch cd] %将４个系数合并在一起，注意cv 与ch的位置与第四章图4-04区别。 对比ca3, ch3,cv3, cd3与AC3，可知矩阵算法正确。细心的同学也许能发现我们的矩阵位置和第４章图４-０４(a)本文图1解释的不同，这是因为dwt2采用的算法与我们不同，但这并不影响小波分解系数的使用。详细解释见附录1。 图1、多级分解的各个系数矩阵 其重构过程只需要两步：（参照第三章24页的公式，但是我们只采用单级重构） AR1=M1*AC1; A=AR1*M1' ; 问题3 如何进行图像的1级非标准分解和重构？ 如果是进行３级非标准分解，那么只需要对前一次系数矩阵的的左上角的近似系数图-04（a）中的a进行，如对是３级分解： M1=haarmatrix(8,8,1) ; AR1=A*M1; 　 %这就是第1次逐行分解的图像 AC1=M1'*AR1; 　 %这就是第１次逐列分解的图像， M2=haarmatrix(4,4,１); AR2=AC1(1:4,1:4)*M2; %第２级行分解，只对左上角低频近似系数部分 AC2=M2'*AR2; 　　　%第２级列分解 M3=haarmatrix(2,2,１); AR3=AC2(1:2,1:2)*M3; %第3级行分解，只对左上角低频近似系数部分 AC3=M3'*AR3; 　　 %第２级列分解 执行MATLAB提供的２维多级小波分解函数wavedec2 [C,L]=wavedec2(A,3,'haar'); ca3=appcoef2(C,L,'haar',３); [ch3,cv3,cd3] = detcoef2('all',C,L,3) AC3相当于上图中(b)左上角的最小的[a3 v3 h3 d3]，AC2是[c2 v2 h2 d2]，AC1是[c1 c2 c3 c4]。那么如果合成图4-04(c)这样的图像只需要把AC3替代AC2的a2 ,在把AC2替代AC1的a1就可以，相应的代码： AC2(1:2,1:2)=AC3; AC1(1:4,1:4)=AC2; 如何由3级分解系数矩阵重构出原始矩阵？参照1级非标准重构,注意此时从第3级向第1级重构，即首先由A3中的a h v d四个系数重构出A2的低频部分a,然后重复这种方法： AR3=M3*AC3; A3=AR3*M3'； AC2(1:2,1:2)=A3； %重构的A3就是A2的低频近似系数，如图4-04(b)的[a h v d] AR2=M2*AC2; A2=AR2*M2'； AC1(1:4,1:4)=A2； %重构的A2就是A1的低频近似系数，如图4-04(a)的a AR1=M1*AC1; A1=AR1*M1'； %A1就是A 注意：上述分解与重构过程与老师第3章中的分解重构过程在一级分解的情况下相同，在第2级分解时就发生了变化。上述分解过程使用了与MATLAB的wavedec2近似的多分辨率分析(multi-resolution analysis), 产生的矩阵系数除了h与v的位置对调外，二者结果完全相同。然而采用第3章21-24页介绍的方法（即使是使用23页的规范化系数矩阵）不能生成与wavedec2相同的矩阵，正是因为如此，采用设置阈值来压缩图像时会与使用wavedec2的结果不同。原因分析及对应的24页的算法的演示程序请参见附录2。（请林老师务必核对）。

suffer · 发表于 2005-8-11 20:36

三种实现非标准分解的方法： 第一种：2.2.3介绍的矩阵A变成 256X256的照片矩阵， 1:2的矩阵下标变成1:64, 1:4 变成1:128, 其它不变，保存系数 AR1,AC1,AR2,AC2,AR3,AC3为PNG格式，就是图3-26所需要的照片，只不过下1级没有显示上一级的细节系数（h v d）部分。使用合成技术类似这样的语句c=[ca cv ; ch cd]应该很方便得出。当然采用附录2中的程序可以更直接，更快速地实现老师的范例。 第二种：按照第3章伪代码形式使用MATLAB一维离散小波变换函数dwt编写。 dwt的用法: [ca cd] = dwt(I, wname); I是你要分解的一维向量，wname是字符串，使用单引号扩起来，如’db9’或’haar’。返回系数ca是低频近似系数，cd是高频细节系数。 由于db9的滤波系数长度为18，dwt需要扩展以减少边界误差，重构函数idwt能够保证恢复到原来的尺寸。所以可以采用不同的扩展边界模式来使用，缺省值是’sym’,它在分解后会使ca和cd的长度变为 (LI + LF-1)/2, LI是向量I的长度，LF是小波滤波系数的长度。也可以改变边界扩展模式，如使用’per’,具体方法是[ca cd]=dwt(I,wname,’mode’,’per’).它的ca 和cd的长度为 LI/2,但是它可能带来边界误差。 我编写了2维非标准分解nsdwt2和重构nsidwt2函数，可参见附录3。 第三种：直接使用MATLAB卷积函数conv2，即仿照dwt2的函数编写。 此处不再详细讨论，有兴趣者参见附录1 这三种方式在MATLAB下计算速度最快的是第一种，最慢的是第二种，第三种速度与第一种速度比相差不大。在笔者PIII733上运行分解wbarb图像过程所用时间（使用TIC,TOC）分别为 0.10, 0.28, 7.23秒。 常见问题： 如何读取、转换和保存图像文件？ [X, map] = imread(‘文件名’) ; %文件名必须带后缀。 使用size(X)来判断图像矩阵的维数，如果是256 256就是2维的，那么map就是它的色板文件，可以使用imshow(X,map)来显示它。如果256 256 3就是3维真彩色图像（True color Image），使用imshow(X)来显示。一般我们把2维的彩色图片称为伪彩色图像（Indexed color image）。 一般来说，不管图像文件是不是PNG格式，最好首先使用MATLAB的imwrite方式重新保存： 对于真彩色文件，imwrite（X, ‘新文件名’） ; 对于伪彩色文件，imwrite（X, map,‘新文件名’）; 由imread读出的图像文件采用uint8的类型，必须转换成double类型之后方可进行运算。转换方法： X=doubl(X); 经过小波变换和重构的X是double类型，在保存图像前需要转换成uint8类型，如果处理不好，会产生截断误差，许多同学反映阈值为0时的文件大小和原始文件大小不一致可能就是出现在这里，转换方法： Y=uint8(round(X)); 这里的round是就近取整，这样7.99999999就会变成8而不会被截断变为7。 如何处理Indexed color 文件？ 原则说，只要可以转换成double类型我们就可以开始小波变换。但是由于伪彩色图片的色板map保存了每个象素值对应的r g b的颜色值，因此，如果map中的颜色值不连续，就不能保证小波转换的效果。可以使用imshow(X,map);colorbar;来显示图片的色板是否连续。如果不连续，可以使用MATLAB提供的方法转换成灰度图像处理。这是林老师在“教师答疑”中贴出的MATLAB帮助文件Image Conversion，我把它转换成一个小程序，可以修改使用。Anna指出ind2gray和ind2rgb也能完成同样的任务。 [X,map]=imread('你的伪彩色图片文件'); X=double(X); if min(min(X))==0 X=X+1; end image(X) title('Original Color Indexed Image') colormap(map); colorbar R = map(X,1); R = reshape(R,size(X)); G = map(X,2); G = reshape(G,size(X)); B = map(X,3); B = reshape(B,size(X)); Xrgb = 0.2990*R + 0.5870*G + 0.1140*B; n = 255; % Number of shades in new indexed image X = round(Xrgb*(n-1)) + 1; map2 = gray(n); image(X), title('Processed Gray Scale Indexed Image') colormap(map2), colorbar baboon= X; map = map2; imwrite(X,map,'gray.png'); %保存为灰度图像，文件名为gray.png 如何处理真彩色文件？ 有几种方法，在MATLAB的Image Conversion中有叙述，但许多同学采用了单独处理它的3个颜色的2维矩阵。设X为三维矩阵(256,256,3) , X(:,:,1)代表红颜色的2维矩阵 X(:,:,2)代表绿2维矩阵, X(:,:,3)代表绿2维矩阵。 [X, map]=imread(‘真彩色文件’); r=double(X(:,:,1)); %r是256 x 256的红色信息矩阵 g=double(X(:,:,2)); %g是256 x 256的绿色信息矩阵 b=double(X(:,:,3)); %b是256 x 256的兰色信息矩阵 %开始对r g b分别作小波变换，并分别形成各自的小波系数矩阵 …… %把3个变换的系数矩阵合并成图像文件 Y(:,:,1)=uint8(round( r )); Y(:,:,2)=uint8(round( g )); Y(:,:,1)=uint8(round( b )); 可以参照3.1的simplecmp.m程序来使用。 如何完成使用小波变换压缩图像的任务 由于压缩的目的在于比较使用不同小波压缩后重构图像的失真度视觉效果和使用PNG文件保存时的文件大小，如果我们自己编写的小波变换程序存在小的遗漏，可能对压缩结果判断错误，所以我认为至少应当使用MATLAB为我们提供的压缩函数记录结果，以便与自己设计的阈值处理程序进行比较。一般来说，两种处理结果应该差别很小，甚至无差别。 直接使用wdencmp: 这是我使用wdencmp编写的简单压缩处理程序simplecmp，可以直接再运行时输出图像、0的个数、PNG文件的大小以及计算时间。如果你的图片为 a.png,不管是真彩色、伪彩色，使用： simplecmp(‘a.png’,’haar’,3) %进行haar小波的三级分解压缩合成 simplecmp(‘a.png’,’db9’,3) %进行db9小波的三级分解压缩合成 下面是具体程序清单。我只对高频细节系数部分进行硬阈值设置。 function simplecmp(fname,wname,level) [rgb,map]=imread(fname); fig=figure;colorbar;axis on;axis equal; set(fig,'position',[10 20 790 580],'name','压缩程序演示') if length(size(rgb))==3 rgbcmp(rgb,wname,level); else indcmp(rgb,map,wname,level); end %--------------------------------------------------- function rgbcmp(rgb,wname,level) %这是压缩真彩色图像 rgb=double(rgb); THR=[0 5 10 20]; otxt=sprintf('%s_zeros.txt',wname); fid=fopen(otxt,'w'); fprintf(fid,'%20s%15s%15s%15s\n','文件名','大小','阈值','0数'); for T=THR TIC; [rgb(:,:,1),cxc,lxc,perf0,perfl2]=wdencmp('gbl',rgb(:,:,1),wname,level,T,'h',1); num0=length(find(abs(cxc) < 0.0000001)); [rgb(:,:,2),cxc,lxc,perf0,perfl2]=wdencmp('gbl',rgb(:,:,2),wname,level,T,'h',1); num0=num0+length(find(abs(cxc) < 0.0000001)); [rgb(:,:,3),cxc,lxc,perf0,perfl2]=wdencmp('gbl',rgb(:,:,3),wname,level,T,'h',1); num0=num0+length(find(abs(cxc) < 0.0000001)); x=uint8(round(rgb)); oname=sprintf('%s_%d.png',wname,T); imwrite(x,oname); tmp=imfinfo(oname); fs=tmp.FileSize; fprintf(fid,'%20s%15d%15d%15d\n',oname,fs,T,num0); e_t=TOC; sTitle=sprintf('%s阈值%d的文件大小%d,0数%u,用时%f,任意键继续...',wname,T,fs,num0,e_t); image(x) title(sTitle); pause end fclose(fid); %--------------------------------------------------- function indcmp(x,map,wname,level) %这是压缩伪彩色图像 THR=[0 5 10 20]; x=double(x); otxt=sprintf('%s_zeros.txt',wname); fid=fopen(otxt,'w'); fprintf(fid,'%20s%15s%15s%15s\n','文件名','大小','阈值','0数'); for T=THR TIC; [y,cxc,lxc,perf0,perfl2]=wdencmp('gbl',x,wname,level,T,'h',1); num0=length(find(abs(cxc) < 0.0000001)); y=uint8(round(y)); oname=sprintf('%s_%d.png',wname,T); imwrite(y,map,oname); tmp=imfinfo(oname); fs=tmp.FileSize; fprintf(fid,'%20s%15d%15d%15d\n',oname,fs,T,num0); e_t=TOC; sTitle=sprintf('%s阈值%d的文件大小%d,0数%u,用时%f,任意键继续...',wname,T,fs,num0,e_t); imshow(y,map) title(sTitle); pause end fclose(fid); 自己编写阈值设置函数： 可以使用wthresh(X,’h’,T)，X是要进行阈值设置的矩阵，’h’表示使用硬阈值处理方式，阈值T=5，如对X=[6,-6 5 2 -2]进行阈值设置的结果是[6 6 0 0 0]，如果设置软’s’，结果是 [1 1 0 0 0],大于T的，也被缩小。 注意:低频系数是对图像重构质量最重要的系数，一般不需要设置。

suffer · 发表于 2005-8-11 20:36

仍需探讨的问题：为什么使用PNG存储经小波变换后的重构图像变大？ 我曾在清华大学的多媒体课程的教师答疑中写了“老师：尊重事实:DB9阈值10的PNG文件就是比原文件大”和“续一：尊重事实:DB9阈值10的PNG文件就是比原文件大”，在林老师的鼓励和指导下，我进行了继续试验、分析，与刘赵璧（Anna）同学进行了探讨，并得到了Lily（姓名还不知道）同学的帮助，同时同学们也做了各自不同的实验，现在的实验结果可以说基本上比较明确，那就是有些图像就是会变大，这与图像的种类、纹理等密切相关。林老师曾经鼓励我去研究一下PNG的压缩方法，无奈我资质不够，至今在这方面的进展不大。由于临近期末考试，作业也要抓紧，所以我暂且将没有搞明白的内容搁置，待寒假期间再进行，希望对这些问题有各种看法也有兴趣研究的同学对此发表意见。以下是我最近试验、分析和阅读到的一些相关信息。 试验结果 我首先根据老师第三章的Haar矩阵算法推演出DB9的系数矩阵，并实现了分解重构及阈值处理程序，对几种照片进行了比较，然后使用3.1节的simplecmp进行了相同照片的实验，结果相当一致。细小差别是因为我的程序对边界的扩展与MATLAB不一样，在设置阈值后引起了边界上小部分不一致造成的。 表一：真彩色图像百合花的处理结果 阈值 Png Haar(Mat/Mine) 0数 Haat (Mat/Mine) PNG Db9 (MAT / Mine) 0数 Db9 (MAT / Mine) 95973 / 95973 0 95973 / 95973 27524 / 24268 95973 / 95973 27 / 9 5 74552 / 74292 135838 / 136063 101882 / 101992 167412 /165662 10 51976 / 51504 163423 / 163741 98411 / 98861 199200 / 195730 20 32474 / 32346 180167 / 180267 92295 / 93660 220629 / 217214 从对比表中我们能够看到2个程序的结果相当一致，因此，我不再给出两种程序的对比，而是使用simplecmp直接处理的结果说明。 将百合花图像使用[I,map]=rgb2ind(x,255);转换成为彩色图像处理，在将伪彩色图像转换为连续变换的灰度图像（如2.4常见问题中讨论的方法）进行处理： 表二：百合花的伪彩色图像和处理后的灰度（gray）图像的处理结果 阈值 Png Haar(Index/Gray) 0数 Haar(Index/Gray) PN Db9(Index/Gray) 0数 Db9 (Index/Gray) 48535 / 43235 0 48535 / 43235 6096/ 7430 48535/ 43235 18 / 22 5 53207 / 36450 9473/ 43626 60362 / 49927 7009 / 52852 10 58025 / 23602 13362/ 54344 64916 / 47813 13202 / 65881 20 60193/ 14347 21948 / 60039 66020/ 46014 24468 / 73494 其他伪彩色与进行加工的灰度图的结果与此完全一致，这也就说明了如果伪彩色文件的色板不是单调性递增就不适合小波分解。“The color bar to the right of the image is not smooth and does not monotonically progress from dark to light. This type of indexed image is not suitable for direct wavelet decomposition with the toolbox and needs to be preprocessed.”。我对 Facets进行同样的实验，结果与此一致。这种处理的结果可以从图像象素值的连续性来理解。这是处理与不处理的图像的中间一行的数据图。另外，不连续的图像质量在压缩后会被极大地破坏 图2伪彩色文件变化前后的第128行数据的连续性情况对比 分析 多种试验图片基本能够反映类似的结果，虽然Indexed Color image有时令Haar小波的分解重构图像出现增大现象，单经过处理之后，这种现象就会消失。然而对于DB9可以看到无论真彩色还是处理后的灰度图像都在阈值5 10处超过原始图像的大小，能不能因此得出DB9不适合进行图像压缩的结论呢？有一些同学确实这样认为，但我认为这种观点因为忽略了如何利用小波进行压缩和还原的过程，这也正是第四章老师为我们讲述的那些编码算法而造成的。在推荐材料[1]中也有类似的说明。 图3、JPEG2000的基本结构 看一下上图就可以明白为什么PNG不能衡量小波压缩的效率问题。上图的图像原始数据首先经过正变换（Forward Transform）就是小波变换的得到小波系数，变换的小波系数经过阈值处理后进行量化，编码后得到压缩的图像文件JPEG2000，如果你没有JPEG2000的显示程序，那么你就不能看到它。它的显示程序就是由解码器从压缩数据中解出编码，进行反量化，得到小波系数，再实施逆变换（Inverse Transform）就是小波系数重构。最终得到图像的原始数据。 因此衡量小波变换的效率是应该看你选择的小波能不能分解出适合“编码器”压缩的小波系数，这种编码器不是PNG的LZ77，因为LZ77压缩小波分解系数的效率不是最好的。这种高效编码器在第四章可以找到。 那么我们存储 PNG文件的目的是什么呢？我认为压缩与去噪(de-noising)是同一种方法的两种提法。他们都使用了设置阈值的方法。我们可以仔细分析经过重构的PNG图片的质量来体会这种消除噪音的效果，也可以评定小波压缩后的图片的视觉质量，同时PNG的文件大小也可以让我们从LZ77算法的本质来理解小波变换压缩后的重构图像的内容变化情况。比如，我们可以从表2中的灰度图像在haar变换取阈值20时出现块状象素，文件大小变为14347，而db9却为46014，超过原始的PNG大小，但并不出现块状而是具有波状的特征。这本身说明了采用Harr小波压缩或去噪后重构的图像中相同的‘串’增多，便于PNG方式压缩，而db9则在相同阈值的情况下不会象Haar那样制造‘马赛克’，说明了它的平滑性，这也能帮助我们理解小波的特性。当然，当阈值继续增加后，超过某一界限，即使DB9也仍然会使PNG文件大小减小。这本身也就是由双尺度（Dyadic）小波变换的两种滤波器决定的。低通滤波结果相当于平均值，高通滤波结果相当于差值，差值能够保证重构图像的细节部分丢失最小，如果差值部分被阈值略去的过多，细节就会越来越少，平均意义的值约来越多，直到多到某一个临界值时（该图像的阈值取到40），重构的图像也可能出现较多的相同数字串，这就会提高PNG的压缩结果。下图是我对Haar（蓝色）小波取阈值为20，db9（红色）小波阈值取40时第128行1：32列的数据曲线与原始数据（黑色）曲线的对比。可以看出也db9在阈值=40时出现了较多的平均值，但比haar在阈值=20时的曲线要少的多。 图4、haar（蓝色）和db9（红色）压缩后重构图像的第128行，1：32列的数据曲线 不过，MATLAB给我们提供了量化的方法来决定如何选取阈值。在HELP Wavelet Toolbox : Advance Concepts: Choosing the Optimal Decomposition中提到了几种利用“熵”的概念来衡量如何选取合适的分解级。感兴趣的同学还可以参看wentropy, wdcbm2, wpdec的帮助。文献[1]中也提到了衡量压缩质量的客观化方法MSE,PNSR并指出小波的重构滤波器的长度越长，形状越规则越能够提供良好的压缩性能。 上面对PNG的讨论因为没有足够的算法分析和程序解读，同时也没有准确的试验数据，因此只能作为猜测。但衡量小波压缩效率的方法我坚持认为不能以PNG文件大小来解说，如果采用图像文件大小来衡量，应该以JPEG2000来衡量。 　 致谢： 感谢林老师的多次指导，感谢广东中山站的刘赵璧同学的帮助和建议并校对修改了许多错误，Lily（我还不知道她的姓名）为我提供的宝贵资料，以及曹红军、李澄钰、郑南、施吉鸿等在网络上共享他们程序的同学。感谢张子亮以及网上许多同学的讨论对我的启发。 　 　 附录： 1、dwt2的分解系数与第3章使用矩阵方法的分解系数不一致的原因 第3 章小波与小波变换 procedure NonstandardDecomposition(C: array [. . . , . . . ] h h 1 1 of reals) C=C/h while h > 1 do for row 1 to h do DecompositionStep(C [row, 1 . . . h]) end for % C is changed to Average and Difference % A-Left half, D-right half for col 1 to h do DecompositionStep(C [1 . . . h, col]) end for % C is changed to ca ch,cv,cd, ca-Top Left, % ch-bottom left,cv-top right, cd --bottom right h=h/2 end while end procedure 1）两个DecompositionStep都处理(C[])，那么C[]在row处理后，原来的C已经改变为行处理后的。应当声明。 2）我仔细阅读了教材，并使用MATLAB编写上述算法，发现老师没有说明产生了几个系数。这是我分析MATLAB的dwt,dwt2,idwt2之后并亲自编写自己的nsdwt2和nsidwt2时发现的。您的原文“首先对图像的每一行计算像素对的均值和差值，然后对每一列计算像素对的均值和差值。”设平均产生结果为A,差为D,试想每一行的平均值是h/2个A,差值是h/2个D,那么仍然存储在C中（两者都假设下采样），此时仍然使用DecompositionStep(C [1 . . . h, col]),就会出现h/4个AA,h/4个AD(A的均值与差值）,h/4个DA,h/4个DD（A的均值与差值），这实际上就是说，2维与1维只有A与D两个系数相比，产生了4个系数，这也就是使用[ca ch cv cd]=dwt2(C,wname)的得到的低频系数，水平、垂直、对角系数。对于最后生成的矩阵C,左上ca,左下ch,右上cv,右下cd。虽然您在前面章节介绍了左上角矩阵AA为低频，其余均为高频细节系数，但如果使用dwt编写上述算法时，很容易只考虑低频部分而丢掉cv,cd。 ------------------ dwt2源程序: 前面的源代码略去…… 首先扩展边界，sym-->sizeEXT=(18-1), per-->18/2 y = wextend('2D',dwtEXTM,x,[sizeEXT,sizeEXT]); 在进行低频Lo_D系数的每行的过滤得到z，这相当于行变换，只不过没有下采样，就是上面的A部分 z = wconv('row',y,Lo_D); 然后再用行变换得低频矩阵z再一次进行列方向低频滤波卷积下采样，得到的低频部分就是 ca=a,对z得高频滤波就是 ch=h,再对比上面对AA,AD a = convdown(z,Lo_D,sizeKEPT,shift); h = convdown(z,Hi_D,sizeKEPT,shift); 在进行高频Hi_D系数的每行的过滤得到z，这相当于行变换，只不过没有下采样，就是上面的D部分 z = wconv('row',y,Hi_D); 然后再用行变换得低频矩阵z再一次进行列方向Lo_D滤波卷积下采样，得到的高频中的低频部分就是 cv=v,对z得高频滤波就是cd=d对角部分,再对比上面对DA,DD v = convdown(z,Lo_D,sizeKEPT,shift); d = convdown(z,Hi_D,sizeKEPT,shift); 如果你想使用dwt2确得到图3-26那样的图像，比使用dwt要快的dwt2中的关键代码的解释： type dwt2.m 与下面的源代码解释相对照，稍加修改，然后换一个函数名称，作业就完成了。 可行的方法：对第一个z进行下采样 a1=dyaddown(wkeep(z,sizeKEPT))就是图3-26的第一次行分解的低频部分，对第2个的z进行上述过程得到d1就是图3-26的高频部分。然后合成图像： img3_26=uint8(round([a1 d1]); 这个速度绝对与dwt2一致，因为它就是dwt2. 　 2、采用第3章21-24页介绍的方法不能生成与wavedec2相同的矩阵。 原来认为老师第三章的haar的2级，3级小波矩阵就只对低频（ca)部分分解,实际上不是，ch部分也参加了计算。但是因为矩阵是正交的，重构没有问题，但压缩可能就会出现问题。举例：4X4 A 与 2级系数 M A=[ a11 a12 a13 a14 M=[ 1/2 1/2 0 0 a21 a22 a23 a24 1/2 -1/2 0 0 a31 a32 a33 a34 0 0 1 0 a41 a42 a43 a44 ] 0 0 0 1 ] 此时a11 a12 a21 a22为ca; a31 a32 a41 a42是ch; A*M= [a11*1/2+a12*1/2 a11*1/2-a12*1/2 a13 a14 a21*1/2+a22*1/2 a21*1/2-a22*1/2 a23 a24 ------------------------------------------------- a31*1/2+a32*1/2 a31*1/2-a31*1/2 a33 a34 a41*1/2+a42*1/2 a41*1/2-a41*1/2 a43 a44

suffer · 发表于 2005-8-11 20:36

由此，ch参加了运算。下面是根据老师第3章21-24页编写的标准和非标准分解重构程序。可以看出它分解的高频部分越来越亮0的数目不是越来越多。但如果将X的赋值语句的注释符号%去掉，可以看到与老师3.5节一样的结果。 function showhaardecrec % This program is completely according Professor Linfuzong's additional materials % chapter 3,section 3.5 2-D haar wavelet transform and synthesis. % Here we use the .mat file wbarb which should be in the directory wavedemo % It use function makehaarmatrix make the haar wavelet transform matrix at 3 different level. % and then 1)showing standard decomposition/reconstruction process and display the result % in a figure 2) showing non-standard decomposition process and display the result in a % figure 3) showing direct wavelet transform/synthesis using the product matrix of 3-levels % matrix and display the image. load wbarb; colormap(map); %X = [ 64 2 3 61 60 6 7 57 %9 55 54 12 13 51 50 16 %17 47 46 20 21 43 42 24 %40 26 27 37 36 30 31 33 %32 34 35 29 28 38 39 25 %41 23 22 44 45 19 18 48 %49 15 14 52 53 11 10 56 %8 58 59 5 4 62 63 1 ] X=double(X); maxX=max(max(X)); X=wcodemat(X,maxX+1); fig=figure; set(fig,'position',[0 0 800 600]); image(X); TIC [rows,cols]=size(X); M1=makehaarmatrix(rows,cols,1,1); M2=makehaarmatrix(rows,cols,2,1); M3=makehaarmatrix(rows,cols,3,1); S1=X*M1; %standard dec row S2=S1*M2; %standard dec row S3=S2*M3; %standard dec row S4=M1'*S3; %standard dec col S5=M2'*S4; %standard dec col S6=M3'*S5; %standard dec col uiwait(msgbox(sprintf('完成标准分解的计算时间 %f',TOC))); subplot(2,3,1);image(wcodemat(S1,maxX+1));title('standard row dec 1'); subplot(2,3,2);image(wcodemat(S2,maxX+1));title('standard row dec 2'); subplot(2,3,3);image(wcodemat(S3,maxX+1));title('standard row dec 3'); subplot(2,3,4);image(wcodemat(S4,maxX+1));title('standard col dec 1'); subplot(2,3,5);image(wcodemat(S5,maxX+1));title('standard col dec 2'); subplot(2,3,6);image(wcodemat(S6,maxX+1));title('standard col dec 3'); TIC S1=X*M1 ; %nsd row S2=M1'*S1; %nsd col S3=S2*M2 ; %nsd row S4=M2'*S3; %nsd col S5=S4*M3 ; %nsd row S6=M3'*S5; %nsd col uiwait(msgbox(sprintf('完成非标准分解的计算时间 %f',TOC))); subplot(2,3,1);image(wcodemat(S1,maxX+1));title('non-standard row dec 1'); subplot(2,3,2);image(wcodemat(S2,maxX+1));title('non-standard col dec 1'); subplot(2,3,3);image(wcodemat(S3,maxX+1));title('non-standard row dec 2'); subplot(2,3,4);image(wcodemat(S4,maxX+1));title('non-standard col dec 2'); subplot(2,3,5);image(wcodemat(S5,maxX+1));title('non-standard row dec 3'); subplot(2,3,6);image(wcodemat(S6,maxX+1));title('non-standard col dec 3'); TIC W=M1*M2*M3; %T=((X*W)'*W)'; T=W'*X*W; %T=(M3'*(M2'*(M1'*X*M1)*M2)*M3) ; rA=wcodemat(((W')^(-1)*T*(W)^-1),maxX+1); %rA=((M1^(-1))'*((M2^(-1))'*((M3^(-1))'*T*(M3)^(-1))*(M2)^(-1))*(M1)^(-1)); uiwait(msgbox(sprintf('完成直接矩阵分解重构的计算时间 %f',TOC))); subplot(2,1,1);image(T);title('using multiplied 3-level matrix dec directly'); subplot(2,1,2);image(rA);title('using multiplied 3-level matrix rec directly'); %-------------------------------------------------------------------------- function [M]=makehaarmatrix(rows,cols,level,flag) M=zeros(rows,cols); if flag== 0 sv=0.5; %non-normalized else sv=1.0/sqrt(2.0); % normalized end half=rows/(2^level) for i= 1:2:half*2 M(i,ceil(i/2)) =sv; M(i+1,ceil(i/2))=sv; M(i,ceil(i/2)+half) =sv; M(i+1,ceil(i/2)+half)=-sv; end if half < rows/2 for i= half*2+1:rows M(i,i)=1; end end M=sparse(M); 　 3、nsdwt2.m与 nsidwt2.m function [ca,ch,cv,cd,car,cdr]=nsdwt2(C,wname) %2-D non standard discrete wavelet transform. %output: ca--coef.of approximation, %ch,cv,cd-horizental,vertical and diagnal coef of detail %car and cdr is coef of appeoximation and detail of C respetively when %do dwt along row direction. %input 2-D matrix to be decomposed and wavelet name to be used [m,n]=size(C); for row = 1:m [car(row,:),cdr(row,:)]=dwt(C(row,:),wname,'mode','per'); end % Notice that after first dwt, length of every row has been changed [m,n]=size(car); for col = 1:n [cal,chl]=dwt(car(:,col)',wname,'mode','per'); ca(:,col)=cal'; ch(:,col)=chl'; [cvl,cdl]=dwt(cdr(:,col)',wname,'mode','per'); cv(:,col)=cvl'; cd(:,col)=cdl'; end %----------------------------------------------- function [X,car,cdr]=nsidwt2(ca,ch,cv,cd,wname) %2-D non standard inverse discrete wavelet transform. %output X reconstruced result, car--when use idwt along col, %car and cdr is coef of appeoximation and detail of C respetively when %do idwt along col direction. %ch,cv,cd-horizental,vertical and diagnal coef of detail %input 2-D matrix to be decomposed and wavelet name to be used [m,n]=size(ca); for col = 1:n cdl=idwt(cv(:,col)',cd(:,col)',wname,'mode','per'); cal=idwt(ca(:,col)',ch(:,col)',wname,'mode','per'); cdr(:,col)=cdl'; car(:,col)=cal'; end [m,n]=size(car); for row=1:m X(row,:)=idwt(car(row,:),cdr(row,:),wname,'mode','per'); end 参考文献： [1] M. Antonini, M. Barlaud, P. Mathieu and I. Daubechies: Image Coding Using the WaveletTransform, IEEE Trans. Image Proc., pp. 205-220, April 1992 小波工具箱函数 函数指令含义 Allnodes 计算树结点 appcoef 提取一维小波变换低频系数 appcoef2 提取二维小波分解低频系数 bestlevt 计算完整最佳小波包树 besttree 计算最佳(优)树 biorfill 双正交样条小波滤波器组 biorwavf 双正交样条小波滤波器 centfrq 求小波中心频率 cgauwavf Complex Gaussian小波 cmorwavf coiflets小波滤波器 cwt 一维连续小波变换 dbaux Daubechies小波滤波器计算 dbwavf Daubechies小波滤波器 dbwavf(W) W='dbN' N=1,2,3,...,50 ddencmp 获取默认值阈值(软或硬)熵标准 depo2ind 将深度-位置结点形式转化成索引结点形式 detcoef 提取一维小波变换高频系数 detcoef2 提取二维小波分解高频系数 disp 显示文本或矩阵 drawtree 画小波包分解树(GUI) dtree 构造DTREE类 dwt 单尺度一维离散小波变换 dwt2 单尺度二维离散小波变换 dwtmode 离散小波变换拓展模式 dyaddown 二元取样 dyadup 二元插值 entrupd 更新小波包的熵值 fbspwavf B样条小波 gauswavf Gaussian小波 get 获取对象属性值 idwt 单尺度一维离散小波逆变换 idwt2 单尺度二维离散小波逆变换 ind2depo 将索引结点形式转化成深度—位置结点形式 intwave 积分小波数

qiuzhuwenti · 发表于 2007-5-14 09:23

向楼主学习

yuansf · 发表于 2008-1-17 16:10

:@) 向楼主学习

yannai · 发表于 2010-3-17 15:30

支持，特别感谢，好好学习和研究！
让我有了学懂小波的强烈动力！

he_sir · 发表于 2010-4-8 22:10

谢谢，学习了。

hanyingsj · 发表于 2010-11-12 10:48

精辟，精辟

hanyingsj · 发表于 2010-11-18 16:22

精辟精辟

账号		自动登录	找回密码
密码			我要加入

[综合讨论] 使用MATLAB设计小波变换程序中的若干问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复：（suffer）使用MATLAB设计小波变换程序中的若干...

回复：（suffer）使用MATLAB设计小波变换程序中的若干...

回复：（suffer）使用MATLAB设计小波变换程序中的若干...

回复：（suffer）使用MATLAB设计小波变换程序中的若干...

请大家帮忙

点评