嵌入式小波零树编码压缩

嵌入式小波零树（EZW）算法的过程详解和Matlab代码（3）解码过程收藏

1、初始化

（1）全局初始化：包括获取扫描次序列表scanorder，初始化符号矩阵flagMat，创建以下几个空表：重要系数重构列表recvalue，量化器编号列表rIlist，上一次解码用到的辅扫描表qu antiflagOld。

（2）每级扫描的初始化：获取本级解码需要的主扫描表scancode和辅扫描表quantiflag。创建解码矩阵DecodeMat，量化符号列表编号qrNum，主扫描表扫描编号scNum。

2、构造逆量化器

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

Q：0 1 1

R：1 1 1

Q：1 1 1 1 0 0 0 1 1 0

R：2 3 3 1 1 1 1 1 1 1

Q：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 0 0 0 1 1 0

R：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 1

Q：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 1……

R：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 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) + quantiflag(r-1)

从而，就可以轻松编写出逆量化器的函数代码：

function [antiQuantiMat,rIlist,quantiflagOld]=antiquantifier(T1,level,rIlist,quantifl ag,quantiflagOld)

antiQuantiMat=[];

maxInterValue=2*T1;

threshold=T1/2^(level-1);

intervalNum=maxInterValue/threshold-1;

for i=1:intervalNum

antiQuantiMat=[antiQuantiMat;threshold*(i+0.25) threshold*(i+0.75)];

end

rIlen=length(rIlist);

flaglen=length(quantiflag);

for r=1:rIlen

rIlist(r)=2*rIlist(r)+quantiflagOld(r);

end

for f=rIlen+1:flaglen

rIlist(f)=1;

end

quantiflagOld=quantiflag;

3、提高上一级解码的重要系数的重构精度

这里主要是根据逆量化器生成的重构值矩阵antiQuantiMat和逆量化器编号列表rIlist来更新重要系数列表recvalue的值，并将更新数据存入解码矩阵DecodeMat中。这四个矩阵（列表）的每一行（个）元素都是一一对应的关系，更新就很容易实现，代码如下：

function [DecodeMat,recvalue,qrNum]=updateRecvalue(DecodeMat,recvalue,qrN um,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=-qValue;

end

recvalue(i,1)=qValue;

DecodeMat(recvalue(i,2),recvalue(i,3))=qValue;

qrNum=qrNum+1; % recvalue矩阵的第2、3列储存的是对应于DecodeMa t中的行、列号（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=qrNum+1;

end

end

end

4、解码

解码过程其实与主扫描过程非常类似，按照扫描次序scanorder扫描解码矩阵DecodeMat：（1）首先检查对应的flagMat是否为“O”或“X”，是则跳过，扫描下一个元素。

（2）然后对照主扫描表scancode，作如下操作：

如果为“P”或“N”，则更新flagMat，根据逆量化器重构DecodeMat相应元素的值，并存入重要系数矩阵recvalue中；

如果为“Z”，则DecodeMat相应元素的值置零；

如果为“T”，则更新flagMat，DecodeMat相应元素的值置零，并使该点所有子孙的flagMat 置为“X”。

（3）本级解码完毕，保存符号矩阵flagMat到解码符号矩阵Decodeflag中，然后使flagMat 中的元素重置为“O”或“Z”，供下一级解码使用。

function [DecodeMat,flagMat,recvalue]=decoding(DecodeMat,flagMat,recvalue,an tiQuantiMat,quantiflag,rIlist,scanorder,scancode,scNum,qrNum)

global row col

for r=1:row*col

if flagMat(scanorder(r,1),scanorder(r,2))=='O'

continue;

elseif flagMat(scanorder(r,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=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 scancode(scNum)=='N'

flagMat(scanorder(r,1),scanorder(r,2))='N';

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=-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 scancode(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(rc,2))~='O'

flagMat(chTree(rc,1),chTree(rc,2))='X';

end

end

scNum=scNum+1;

end

end

end

Decodeflag=flagMat;

for i=1:row

for j=1:col

switch flagMat(i,j)

case {'P','N'}

flagMat(i,j)='O';

case {'T','X'}

flagMat(i,j)='Z';

end

end

end

嵌入式小波零树（EZW）算法的过程详解和Matlab代码（4）运行示例收藏

最后，我们以一幅64*64的Lena局部图像来演示EZW算法的编解码效果。

首先给出相关的编解码主程序。

function y=ezw(Im,imDim,codeDim,DecodeDim)

global row col

[row,col]=size(Im);

% ----- Wavedec -----%

DecIm=mywavedec2(Im,imDim);

% ----- Select threshold -----%

MaxDecIm=max(max(abs(DecIm)));

T=zeros(1,codeDim);

T(1)=2^flor(log2(MaxDecIm));

for i=2:codeDim

T(i)=T(i-1)/2;

end

% ----- EZW coding -----%

[CodeList,LenSubCL,QuantiFlagList,LenSubQFL]=ezwcode(DecIm,T,codeDim);

% ----- EZW decoding -----%

DecodeMat=ezwdecode(DecIm,T(1),DecodeDim,CodeList,LenSubCL,QuantiFlagLis t,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 col

scanlist=morton(Mat);

flaglist(1:row,1:col)='Z';

imptvalue=[];

imptflag=[];

% ----- Intializing EZW coding output variables -----%

CodeList=[];

LenSubCL=[];

QuantiFlagList=[];

LenSubQFL=[];

% ----- Coding loop -----%

for d=1:codedim

[imptvalue,imptflag,scancode,scanflag,flaglist]=mainscan(Mat,scanlist,flaglist,i

mptvalue,imptflag,threshold(d));

[quantilist,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)];

end

function DecodeMat=ezwdecode(Mat,T1,decodeDim,CodeList,LenSubCL,QuantiFla gList,LenSubQFL)

global row col

recvalue=[];

rIlist=[];

quantiflagOld=[];

scanorder=listorder(row,col,1,1);

flagMat(1:row,1:col)='Z';

for level=1:decodeDim

scancode=CodeList(1:LenSubCL(level));

CodeList=CodeList(LenSubCL(level)+1:end);

quantiflag=QuantiFlagList(1:LenSubQFL(level));

QuantiFlagList=QuantiFlagList(LenSubQFL(level)+1:end);

DecodeMat=Mat;

qrNum=1;

scNum=1;

[antiQuantiMat,rIlist,quantiflagOld]=antiquantifier(T1,level,rIlist,quantiflag,qua ntiflagOld);

[DecodeMat,recvalue,qrNum]=updateRecvalue(DecodeMat,recvalue,qrNum,qu antiflag,antiQuantiMat,rIlist);

[DecodeMat,flagMat,recvalue]=decoding(DecodeMat,flagMat,recvalue,antiQua ntiMat,quantiflag,rIlist,scanorder,scancode,scNum,qrNum);

% if level

% button = questdlg('Need the nest decode processing ?','Decode Level',' Yes','No','Yes');

% if button(1)=='N'

% break;

% end

% end

end

（1）图像3级分解、5级编码、5级解码下的效果图。

ezw(Im64,3,5,5);

分解矩阵DecIm的部分数据

40.9531 87.2344 94.2188 161.2031 64.3125 40.4063 37.5781 65.6 875

110.8281 129.3125 91.5625 164.3750 107.1719 50.5469 34.9688 75. 2188

150.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.3906 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 30 82

118 118 94 138 106 30 34 78

102 102 118 134 70 30 38 74

102 86 102 106 30 30 46 78

98 122 134 74 30 34 50 86

运行效果：

（2）图像3级分解、10级编码、8级解码下的效果图。

ezw(Im64,3,10,8);

分解矩阵DecIm的部分数据

40.9531 87.2344 94.2188 161.2031 64.3125 40.4063 37.5781 65.6 875

110.8281 129.3125 91.5625 164.3750 107.1719 50.5469 34.9688 75. 2188

150.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.3906 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的对应数据

40.7500 87.2500 94.2500 161.2500 64.2500 40.2500 37.7500 65.7 500

110.7500 129.2500 91.7500 164.2500 107.2500 50.7500 34.7500 75. 2500

150.7500 134.7500 98.2500 164.7500 121.7500 67.7500 28.2500 77. 2500

132.7500 117.7500 91.2500 155.2500 117.2500 53.2500 28.2500 83. 2500

119.2500 118.7500 95.2500 137.2500 105.7500 30.2500 33.7500 79. 7500

100.7500 102.7500 119.2500 132.7500 71.7500 28.7500 38.7500 75. 2500

103.2500 85.2500 102.7500 107.2500 31.2500 28.2500 45.2500 79. 7500

98.2500 120.2500 134.7500 73.2500 29.2500 33.7500 50.2500 86. 2500

运行效果：

EZW算法的过程详解和Matlab代码（1）构建扫描次序表（修正小波树结构）收藏

前段时间，我们讨论了嵌入式小波零树算法的基本原理。（https://www.360docs.net/doc/a1504413.html,/chenyusiy uan/archive/2007/11/16/1888968.aspx）一个多星期过去了，我根据自己对算法的理解，编写出EZW算法的Matlab程序，可以实现图像的任意级别的小波分解和重构、以及任意精度的EZW编解码过程。下面，我们以一幅16*16的Lena局部图像经过3级小波分解后的小波数据为例详细的说明EZW算法的编解码过程，并给出相应的Matlab代码。原始图像如下

分解后的小波图像数据为：

一、在开始编码之前，首先要求出初始阈值T1：

MaxDecIm=max(max(abs(DecIm)));

T=zeros(1,codeDim);

T(1)=2^flor(log2(MaxDecIm));

二、然后是建立小波树结构，构建扫描次序表。这个扫描次序表非常重要，后面的编码、解码过程都要按照扫描次序表逐个处理数据矩阵的各个元素。构建过程如下：用（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）第N~2高频子带(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），给出该点的子孙列表cp 。

function cp=treeMat(r,c) %这个函数是一个递归函数

global row col dim % dim是小波分解级数

HLL=row/2^dim;

WLL=col/2^dim;

if (r<=HLL)&&(c<=WLL)

tp1=[r,c+WLL;r+HLL,c;r+HLL,c+WLL];

cP=[tp1;treeMat(r,c+WLL);treeMat(r+HLL,c);treeMat(r+HLL,c+WLL)];

elseif (r>row/2)||(c>col/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)<=col/2)

t1=treeMat(tp(1,1),tp(1,2));

tm1=[tm1;t1];

t2=treeMat(tp(2,1),tp(2,2));

tm2=[tm2;t2];

t3=treeMat(tp(3,1),tp(3,2));

tm3=[tm3;t3];

t4=treeMat(tp(4,1),tp(4,2));

tm4=[tm4;t4];

end

cP=[tp;tm1;tm2;tm3;tm4];

end

示例，当row=8，col=8，dim=2时，LL-N低频子带的点(1,1) (1,2) (2,1) (2,2) 的子孙分布如下：

由这个小波树列表tree_p，我们可以进一步构建函数childMat()，给出矩阵数据Mat和矩阵任一点的位置（r，c），返回该点的子孙数据列表chMat。

function chMat=childMat(Mat,chRows,chCols)

global row col dim

chPoint=treeMat(chRows,chCols);

chMat=[];

[mRows,mCols]=size(chPoint);

for iRows=1:mRows

chMat=[chMat;chPoint(iRows,1),chPoint(iRows,2),Mat(chPoint(iRows,1),chPoint (iRows,2))];

end

2、构建扫描次序表

本文EZW算法的扫描次序为Morton式，其特征是从（1，1）开始，每4个点组成一个“Z”型扫描单元，从微观到宏观上都是严格的“Z”型结构，可以用递归方法来构建扫描次序表。扫描次序表scanlist由两部分组成，一个是将数据矩阵Mat按照morton扫描次序转换成数据列表matlist，一个是按照扫描次序组成的矩阵各点位置的（r，c）列表lsorder。

function scanlist=morton(Mat)

global row col

matlist=mat2list(Mat);

scanorder=listorder(row,col,1,1);

scanlist=[];

for i=1:row*col

scanlist=[scanlist;i scanorder(i,:) matlist(i)];

end

function mls=mat2list(Mat) % 该函数为递归函数

[r,c]=size(Mat);

if (r==2)&&(c==2)

mls=[Mat(1,1);Mat(1,2);Mat(2,1);Mat(2,2)];

else

M1=Mat(1:r/2,1:c/2);

M2=Mat(1:r/2,c/2+1:c);

M3=Mat(r/2+1:r,1:c/2);

M4=Mat(r/2+1:r,c/2+1:c);

lt1=mat2list(M1);

lt2=mat2list(M2);

lt3=mat2list(M3);

lt4=mat2list(M4);

mls=[lt1;lt2;lt3;lt4];

end

function lsorder=listorder(mr,mc,pr,pc) % 该函数为递归函数

lso=[pr,pc;pr,pc+mc/2;pr+mr/2,pc;pr+mr/2,pc+mc/2];

mr=mr/2;mc=mc/2;

lm1=[];lm2=[];lm3=[];lm4=[];

if (mr>1)&&(mc>1)

ls1=listorder(mr,mc,lso(1,1),lso(1,2));

lm1=[lm1;ls1];

ls2=listorder(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);

嵌入式小波零树（EZW）算法的过程详解和Matlab代码（2）编码过程（0812 27-补充quantifier程序）收藏

首先根据实际需要确定扫描次数D=codeDim，然后进行D次主扫描和辅扫描。

1、初始化

首先获取扫描次序表scanlist，表中每行的数字依次是该元素的（1）扫描序号、（2）对应于数据矩阵中位置的行号r、（3）对应于数据矩阵中位置的列号c、（4）值Mat（r，c），即

scanlist (i) = [ i, r, c, Mat(r, c) ]

编码扫描的对象就是次序表scanlist，而非数据矩阵Mat。然后构建初始符号矩阵flagMat，其中每个元素均为字符“Z”。

2、主扫描

在第d次（1<= d <=codeDim）扫描中，输入阈值为T（d），首先依次对次序表的每一行scanlist(i)按照如下步骤扫描：

（1）如果上一次扫描得到的flagMat中，处于位置（r，c）的符号为“O”，则表示相应的Mat（r，c）是重要的，不需要通过本次扫描再次编码，故作跳过处理。

（2）将scanlist(i,4)（即Mat(r,c)）的绝对值与阈值T(d) 比较，如果是重要的系数，则令flagMat(r,c)的值为P或N，并存入主扫描表scancode中，并将Mat(r,c)的值以及相应的行、列数（r，c）存入重

要数列表imptvalue中，将符号P或N存入重要数符号列表imptfla g中。

（3）如果是不重要的系数，则首先检查flagMat(r,c)，若是符号“X”，则表示这个点是不重要的，作跳过处理不再扫描。

（4）如果点（r，c）处于第一分解级，即r>row/2或c>col/2，没有子孙系数，且其系数是不重要的，则该点flagMat的符号为Z，存入主扫描表scancode中。

（5）如果点（r，c）满足r

至此，扫描结束。然后，将扫描后得到的符号矩阵flagMat转存到sc anflag中，再将flagMat中的字符P、N都置换为O，字符X和T都置换为Z，使得flagMat的元素只有字符O、Z两种，作为输出到下一次扫描中使用。

function [imptvalue,imptflag,scancode,scanflag,flaglist]=mai nscan(Mat,scanlist,flaglist,imptvalue,imptflag,threshold) global row col

scancode=[];

for i=1:row*col