EZW算法的过程详解和Matlab代码

1、精选优质文档-倾情为你奉上精选优质文档-倾情为你奉上专心-专注-专业专心-专注-专业精选优质文档-倾情为你奉上专心-专注-专业 我们以一幅16*16的Lena局部图像经过3级小波分解后的小波数据为例详细的说明EZW算法的编解码过程，并给出相应的Matlab代码。原始图像如下分解后的小波图像数据为：一、在开始编码之前，首先要求出初始阈值T1： MaxDecIm=max(max(abs(DecIm); T=zeros(1,codeDim); T(1)=2flor(log2(MaxDecIm); 二、然后是建立小波树结构，构建扫描次序表。这个扫描次序表非常重要，后面的编码、解码过程都要按照扫描次序表

2、逐个处理数据矩阵的各个元素。构建过程如下：用（r，c）表示数据矩阵上各元素的位置。row、col作为全局变量，表示数据矩阵的行、列数。1、小波树结构的特点：（1）对于LL-N低频子带的点(r,c)，有3个孩子：(r,c+W)、(r+H,c)和(r+H,c+W)，其中W、H分别是LL-N子带的宽和高；（2）第N2高频子带(LH、HL、HH)的点都有4个孩子，即： tp=2*r-1,2*c-1;2*r-1,2*c;2*r,2*c-1;2*r,2*c; （3）第1高频子带的点没有孩子。根据小波树的这个特点，可编写如下“小波树”函数treeMat()，输入矩阵内任一点的位置（r，c），给出该点的子孙列

3、表cp 。function cp=treeMat(r,c) %这个函数是一个递归函数global row col dim % dim是小波分解级数HLL=row/2dim; WLL=col/2dim; if (r=HLL)&(crow/2)|(ccol/2) cP=;else tp=2*r-1,2*c-1;2*r-1,2*c;2*r,2*c-1;2*r,2*c; tm1=;tm2=;tm3=;tm4=; if (tp(4,1)=row/2)&(tp(4,2)1)&(mc1) ls1=listorder(mr,mc,lso(1,1),lso(1,2); lm1=lm1;ls1; ls2=list

4、order(mr,mc,lso(2,1),lso(2,2); lm2=lm2;ls2; ls3=listorder(mr,mc,lso(3,1),lso(3,2); lm3=lm3;ls3; ls4=listorder(mr,mc,lso(4,1),lso(4,2); lm4=lm4;ls4; end lsorder=lso;lm1;lm2;lm3;lm4; len=length(lsorder); lsorder=lsorder(len-mr*mc*4+1:len,1:2); 首先根据实际需要确定扫描次数D=codeDim，然后进行D次主扫描和辅扫描。初始化首先获取扫描次序表scanlist

5、，表中每行的数字依次是该元素的（1）扫描序号、（2）对应于数据矩阵中位置的行号r、（3）对应于数据矩阵中位置的列号c、（4）值Mat（r，c），即 scanlist (i) = i, r, c, Mat(r,c) 编码扫描的对象就是次序表scanlist，而非数据矩阵Mat。然后构建初始符号矩阵flagMat，其中每个元素均为字符“Z”。主扫描在第d次（1= d row/2或ccol/2，没有子孙系数，且其系数是不重要的，则该点flagMat的符号为Z，存入主扫描表scancode中。（5） 如果点（r，c）满足rrow/2且c=threshold if scanlist(i,4)=0 fla

6、glist(scanlist(i,2),scanlist(i,3)=P; scancode=scancode P; imptvalue=imptvalue;abs(scanlist(i,4),scanlist(i,2),scanlist(i,3); imptflag=imptflag P; else flaglist(scanlist(i,2),scanlist(i,3)=N; scancode=scancode N; imptvalue=imptvalue;abs(scanlist(i,4),scanlist(i,2),scanlist(i,3); imptflag=imptflag N;

7、end else if flaglist(scanlist(i,2),scanlist(i,3)=X continue; elseif irow*col/4 scancode=scancode Z; else chImt,chMat=childImportant(Mat,scanlist(i,2),scanlist(i,3),threshold,flaglist); if chImt flaglist(scanlist(i,2),scanlist(i,3)=Z; scancode=scancode Z; else flaglist(scanlist(i,2),scanlist(i,3)=T;

8、scancode=scancode T; rowch,colch=size(chMat); for r=1:rowch if flaglist(chMat(r,1),chMat(r,2)=O flaglist(chMat(r,1),chMat(r,2)=X; end end end end endendscanflag=flaglist; for r=1:row for c=1:col switch flaglist(r,c) case P,N flaglist(r,c)=O; case X,T flaglist(r,c)=Z; end endendfunction chImt,chMat=c

9、hildImportant(Mat,chRows,chCols,threshold,flaglist) global row colchMat=childMat(Mat,chRows,chCols); if max(abs(chMat(:,3)=threshold chImt=1; else chImt=0; end辅扫描辅扫描是对输出符号为P和N的系数进行量化，并将量化符号存入辅扫描表quantiflag中。在主扫描中，将输出符号为P和N的系数的符号信息和数值分别储存在列表imptflag和imptvalue中。辅扫描的具体过程是：构造量化器量化器的构造比较简单，不过关键的问题是要确定每一级

10、扫描中量化器的数目。首先，量化区间的最大区间值应该是初始阈值的2倍，最小值是当前的阈值，故第d级扫描时，总的量化区间为 T(d), 2*T(1) )；而量化间隔为T(d)，故第d级扫描的量化器数目为：Q = 2*T(1) T(d) / T(d) 。确定量化器数目后，就要确定每个量化器中“0”“1”对应的量化值。根据算法原理，第q个量化器中，相应于符号“0”的量化值为：T(d)*(q+0.25)，相应于符号“1”的量化值为：T(d)*(q+0.75)。至此，量化器构造完成。function quantifierMat,threshold=quantifier(T1,level) quantifi

11、erMat=;maxInterValue=2*T1; threshold=T1/2(level-1); intervalNum=maxInterValue/threshold-1; for i=1:intervalNum quantifierMat=quantifierMat;threshold*(i+0.25) threshold*(i+0.75);end辅扫描过程首先要确定这个重要系数属于量化器。根据总量化区间与量化间隔的关系，只要将重要系数除以量化间隔，所得的商就是该重要系数所在的量化器编号。即：rI=flor(imptvalue(j) /threshold);将重要系数减去量化器的起始

12、值，所得的值与量化间隔的二分之一值相比较，小于的话，则量化符号为“0”，大于则量化符号为“1”。然后选择量化重构值，存入重构列表recvalue中。function quantilist,quantiflag,recvalue,quantifierMat=assistscan(imptvalue,dim,T1) quantilist=;quantiflag=;recvalue=; quantifierMat,threshold=quantifier(T1,dim); imRow,imCol=size(imptvalue); for j=1:imRowrI=flor(imptvalue(j)/t

13、hreshold); flag01=imptvalue(j)-rI*threshold; if flag01threshold/2 quantiflag=quantiflag;0; recvalue=recvalue;quantifierMat(rI,1); else quantiflag=quantiflag;1; recvalue=recvalue;quantifierMat(rI,2); endendquantilist=imptvalue(:,1),quantiflag,recvalue,imptvalue(:,2),imptvalue(:,3);编码输出每级扫描完成后，把主扫描输出的

14、主扫描表scanlist存入列表Codelist中，并用一个辅助列表LenSubCL记录这一级主扫描表的表长（在信号传输应用中不一定用得上）；把辅扫描输出的量化符号列表quantiflag存入列表QuantiFlagList中，并用一个辅助列表LenSubQFL记录辅扫描表的表长。以下是经过3次扫描后的编码输出code_level = 1scancode = ZPPPTTTTTTTTTTTTquantiflag = 011code_level = 2scancode = PZTTTZNTTTTTTTZTNTTTZNTTTZZZNZNZZZZNZZZZZquantiflag = code_lev

15、el = 3scancode = NPTTNNZTNTTTNTPNPTTNZNZTTZZTTZZTTNTTTZZZZNZZZPPNZZNZPZPZZZZZZZNZ NZZZNNPZZZZPZZPPZZPZZZNZZZZZZZquantiflag = 1、初始化（1）全局初始化：包括获取扫描次序列表scanorder，初始化符号矩阵flagMat，创建以下几个空表：重要系数重构列表recvalue，量化器编号列表rIlist，上一次解码用到的辅扫描表quantiflagOld。（2）每级扫描的初始化：获取本级解码需要的主扫描表scancode和辅扫描表quantiflag。创建解码矩阵Decod

16、eMat，量化符号列表编号qrNum，主扫描表扫描编号scNum。2、构造逆量化器（1）逆量化器包括两个部分，一是量化值部分，与编码程序中的量化器函数代码相同；二是量化器编号列表rIlist的构造，这个是难点和重点，把它构造出来，才能解码出重要系数的数值。（2）以下是编码过程中经过4次扫描产生的辅扫描表quantiflag和量化器编号列表rIlist：Q：0 1 1R：1 1 1Q：1 1 1 1 0 0 0 1 1 0R：2 3 3 1 1 1 1 1 1 1Q：1 0 1 1 0 1 0 0 1 0 0 1 1 0 0 1 1 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 1

17、 0 0 0 1 1 0R：5 7 7 3 2 2 2 3 3 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1Q：0 1 1 1 0 0 1 0 0 1 1 0 0 0 1 0 0 1 1 1 1 0 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 1 0 1 1 0 0 1 0 1 0 0 1 1 0 0 1R：11 14 15 7 4 5 4 6 7 4 2 3 3 2 2 3 3 2 2 2 3 2 3 2 3 2 2 3 2 3 2 3 3 2 2 2 3 3 2 1 1 1 1 1 1 1 1 1

18、 1 1 1 1 1 1仔细观察可以发现两者之间存在一种函数关系：以第一级扫描产生的三个量化符号（0，1，1）为例，其编号列表为（1，1，1），在随后的第2、3、4级分解后，量化符号相继变化为（1，1，1）、（1，0，1）、（0，1，1），编号列表相继变化为（2，3，3）、（5，7，7）、（11，14，15），这里有一个规律：（2，3，3） = 2*（1，1，1）+（0，1，1）（5，7，7） = 2*（2，3，3）+（1，1，1）（11，14，15）= 2*（5，7，7）+（1，0，1）这样，我们就得到如下公式：设r表示扫描级数，则有： rIlist(r) = 2*rIlist(r-1) +

19、 quantiflag(r-1)从而，就可以轻松编写出逆量化器的函数代码：function antiQuantiMat,rIlist,quantiflagOld=antiquantifier(T1,level,rIlist,quantiflag,quantiflagOld) antiQuantiMat=;maxInterValue=2*T1; threshold=T1/2(level-1); intervalNum=maxInterValue/threshold-1; for i=1:intervalNumantiQuantiMat=antiQuantiMat;threshold*(i+0.2

20、5) threshold*(i+0.75); endrIlen=length(rIlist); flaglen=length(quantiflag); for r=1:rIlen rIlist(r)=2*rIlist(r)+quantiflagOld(r); endfor f=rIlen+1:flaglen rIlist(f)=1; endquantiflagOld=quantiflag;3、提高上一级解码的重要系数的重构精度这里主要是根据逆量化器生成的重构值矩阵antiQuantiMat和逆量化器编号列表rIlist来更新重要系数列表recvalue的值，并将更新数据存入解码矩阵Decode

21、Mat中。这四个矩阵（列表）的每一行（个）元素都是一一对应的关系，更新就很容易实现，代码如下：function DecodeMat,recvalue,qrNum=updateRecvalue(DecodeMat,recvalue,qrNum,quantiflag,antiQuantiMat,rIlist) if isempty(recvalue)rvRow,rvCol=size(recvalue); for i=1:rvRow if quantiflag(qrNum)=1 qValue=antiQuantiMat(rIlist(qrNum),2); if recvalue(i)0 qValue

22、=-qValue; end recvalue(i,1)=qValue; DecodeMat(recvalue(i,2),recvalue(i,3)=qValue; qrNum=qrNum+1; % recvalue矩阵的第2、3列储存的是对应于DecodeMat中的行、列号（r,c） else qValue=antiQuantiMat(rIlist(qrNum),1); if recvalue(i)0 qValue=-qValue; end recvalue(i,1)=qValue; DecodeMat(recvalue(i,2),recvalue(i,3)=qValue; qrNum=qrN

23、um+1; endendend解码解码过程其实与主扫描过程非常类似，按照扫描次序scanorder扫描解码矩阵DecodeMat：（1） 首先检查对应的flagMat是否为“O”或“X”，是则跳过，扫描下一个元素。然后对照主扫描表scancode，作如下操作：如果为“P”或“N”，则更新flagMat，根据逆量化器重构DecodeMat相应元素的值，并存入重要系数矩阵recvalue中；如果为“Z”，则DecodeMat相应元素的值置零；如果为“T”，则更新flagMat，DecodeMat相应元素的值置零，并使该点所有子孙的flagMat置为“X”。（3）本级解码完毕，保存符号矩阵flagM

24、at到解码符号矩阵Decodeflag中，然后使flagMat中的元素重置为“O”或“Z”，供下一级解码使用。function DecodeMat,flagMat,recvalue=decoding(DecodeMat,flagMat,recvalue,antiQuantiMat,quantiflag,rIlist,scanorder,scancode,scNum,qrNum) global row colfor r=1:row*colif flagMat(scanorder(r,1),scanorder(r,2)=O continue; elseif flagMat(scanorder(r,

25、1),scanorder(r,2)=X continue; else if scancode(scNum)=P flagMat(scanorder(r,1),scanorder(r,2)=P; if quantiflag(qrNum)=1 qValue=antiQuantiMat(rIlist(qrNum),2); recvalue=recvalue;qValue,scanorder(r,1),scanorder(r,2); DecodeMat(scanorder(r,1),scanorder(r,2)=qValue; qrNum=qrNum+1; else qValue=antiQuanti

26、Mat(rIlist(qrNum),1); recvalue=recvalue;qValue,scanorder(r,1),scanorder(r,2); DecodeMat(scanorder(r,1),scanorder(r,2)=qValue; qrNum=qrNum+1; end scNum=scNum+1; elseif scancode(scNum)=N flagMat(scanorder(r,1),scanorder(r,2)=N; if quantiflag(qrNum)=1 qValue=-antiQuantiMat(rIlist(qrNum),2); recvalue=re

27、cvalue;qValue,scanorder(r,1),scanorder(r,2); DecodeMat(scanorder(r,1),scanorder(r,2)=qValue; qrNum=qrNum+1; else qValue=-antiQuantiMat(rIlist(qrNum),1); recvalue=recvalue;qValue,scanorder(r,1),scanorder(r,2); DecodeMat(scanorder(r,1),scanorder(r,2)=qValue; qrNum=qrNum+1; end scNum=scNum+1; elseif sc

28、ancode(scNum)=Z DecodeMat(scanorder(r,1),scanorder(r,2)=0; scNum=scNum+1; elseif scancode(scNum)=T flagMat(scanorder(r,1),scanorder(r,2)=T; DecodeMat(scanorder(r,1),scanorder(r,2)=0; chTree=treeMat(scanorder(r,1),scanorder(r,2); rowch,colch=size(chTree); for rc=1:rowch if flagMat(chTree(rc,1),chTree

29、(rc,2)=O flagMat(chTree(rc,1),chTree(rc,2)=X; end end scNum=scNum+1; endendendDecodeflag=flagMat;for i=1:rowfor j=1:col switch flagMat(i,j) case P,N flagMat(i,j)=O; case T,X flagMat(i,j)=Z; endendend 最后，我们以一幅64*64的Lena局部图像来演示EZW算法的编解码效果。首先给出相关的编解码主程序。function y=ezw(Im,imDim,codeDim,DecodeDim) global

30、 row col row,col=size(Im); % - Wavedec -%DecIm=mywavedec2(Im,imDim); % - Select threshold -%MaxDecIm=max(max(abs(DecIm); T=zeros(1,codeDim); T(1)=2flor(log2(MaxDecIm); for i=2:codeDimT(i)=T(i-1)/2; end% - EZW coding -% CodeList,LenSubCL,QuantiFlagList,LenSubQFL=ezwcode(DecIm,T,codeDim); % - EZW deco

31、ding -%DecodeMat=ezwdecode(DecIm,T(1),DecodeDim,CodeList,LenSubCL,QuantiFlagList,LenSubQFL); % - Waverec -%DecIm(1:row/8,1:col/8) DecodeMat(1:row/8,1:col/8) RecIm=mywaverec2(DecodeMat,imDim);function CodeList,LenSubCL,QuantiFlagList,LenSubQFL=ezwcode(Mat,threshold,codedim)global row colscanlist=mort

32、on(Mat); flaglist(1:row,1:col)=Z; imptvalue=;imptflag=;% - Intializing EZW coding output variables -%CodeList=;LenSubCL=;QuantiFlagList=;LenSubQFL=;% - Coding loop -%for d=1:codedimimptvalue,imptflag,scancode,scanflag,flaglist=mainscan(Mat,scanlist,flaglist,imptvalue,imptflag,threshold(d); quantilis

33、t,quantiflag,recvalue,quantifierMat=assistscan(imptvalue,d,threshold(1); % Produce code dataflow CodeList=CodeList,scancode; LenSubCL=LenSubCL,length(scancode); QuantiFlagList=QuantiFlagList,quantiflag; LenSubQFL=LenSubQFL,length(quantiflag); endfunction DecodeMat=ezwdecode(Mat,T1,decodeDim,CodeList

34、,LenSubCL,QuantiFlagList,LenSubQFL) global row colrecvalue=;rIlist=;quantiflagOld=;scanorder=listorder(row,col,1,1); flagMat(1:row,1:col)=Z;for level=1:decodeDimscancode=CodeList(1:LenSubCL(level); CodeList=CodeList(LenSubCL(level)+1:end); quantiflag=QuantiFlagList(1:LenSubQFL(level); QuantiFlagList

35、=QuantiFlagList(LenSubQFL(level)+1:end); DecodeMat=Mat; qrNum=1; scNum=1; antiQuantiMat,rIlist,quantiflagOld=antiquantifier(T1,level,rIlist,quantiflag,quantiflagOld); DecodeMat,recvalue,qrNum=updateRecvalue(DecodeMat,recvalue,qrNum,quantiflag,antiQuantiMat,rIlist); DecodeMat,flagMat,recvalue=decodin

36、g(DecodeMat,flagMat,recvalue,antiQuantiMat,quantiflag,rIlist,scanorder,scancode,scNum,qrNum); % if leveldecodeDim% button = questdlg(Need the nest decode processing ?,Decode Level,Yes,No,Yes); % if button(1)=N% break; % end% endend图像3级分解、5级编码、5级解码下的效果图。ezw(Im64,3,5,5);分解矩阵DecIm的部分数据 40.9531 87.2344

37、94.2188 161.2031 64.3125 40.4063 37.5781 65.6875 110.8281 129.3125 91.5625 164.3750 107.1719 50.5469 34.9688 75.2188150.7813 134.5156 98.0625 164.9688 121.7344 67.5469 28.0938 77.1250 132.9219 117.7656 91.0938 155.3125 117.0781 53.0781 28.1563 83.3750 119.3281 118.5156 95.1875 137.1719 105.8750 30.3

38、906 33.6094 79.6406 100.8906 102.9375 119.2031 132.8438 71.5000 28.6875 38.5313 75.3281 103.2813 85.3750 102.6250 107.4531 31.3594 28.2656 45.3438 79.7813 98.4844 120.0938 134.8438 73.3281 29.1719 33.6094 50.3906 86.0781解码矩阵DecodeMat的对应数据42 86 94 162 66 42 38 66 110 130 90 166 106 50 34 74 150 134 98 166 122 66 30 78 134 118 90 154 118 54