多级树集合分裂SPIHT算法的过程详解与Matlab实现

1、多级树集合分裂(SPIHT)算法的过程详解与Matlab实现上星期我们讨论了EZW算法，很高兴收到了一些朋友的email,对算法进行探讨、交流。这也是我开这个博客的源动力之一，学习就应该开诚布公、交流互助，在探讨中加深对所学知识的理解和掌握。在弄懂了EZW算法原理并用Matlab实现后，我继续学习EZW的改进算法。至今有一周的时间没更新博客、写新文章了，其实就是把时间用在EZW的一个改进算法多级树集合分裂(SetPartitioninginHierarchicalTrees,SPIHT)算法的研究和Matlab实现。由于EZW是SPIHT的基础，所以在EZW算法的Matlab代码的基础上，我很

2、快就完成了SPIHT的代码编写，但最痛苦的是一开始没吃透算法原理，程序在初始分集上出了错，调试了两天找不出根本问题，昨天从头再看一次算法原理，才发现问题所在呵呵，小小的粗心就耽搁了我两三天的时间和精力！问题解决后，就编写程序注释了，上次EZW算法的代码都没写注释，让大家看着辛苦，不好意思哦！好，接下来就开始讨论SPIHT算法的原理，然后给出具体的Matlab代码。一、SPIHT算法与EZW算法EZW算法是一种基于零树的嵌入式图象编码算法，虽然在小波变换系数中，零树是一个比较有效的表示不重要系数的数据结构，但是，在小波系数中还存在这样的树结构，它的树根是重要的，除树根以外的其它结点是不重要的。对

3、这样的系数结构，零树就不是一种很有效的表示方法。A.Said和W.A.Pearlman根据Shapiro零树编码算法(EZW)的基本思想，提出了一种新的且性能更优的实现方法，即基于多级树集合分裂排序(SetPartitioninginHierarchicalTrees,SPIHT)的编码算法。它采用了空间方向树(SOT：spatialorientationtree)、全体子孙集合D(i,j)和非直系子孙集合L(i,j)的概念以更有效地表示上述特征的系数结构，从而提高了编码效率。SPIHT算法能够生成一个嵌入位流(embeddedbitstream),使接收的位流在任意点中断时，都可解压和重构图

4、像，具有良好的渐进传输特性；算法的初始化过程、细化过程类似于EZW算法，它改进了EZW重要图的表示方法，也就是重要系数在表中的排序信息，使得集合的表示更为精简，从而提高了编码效率和图像压缩率。SPIHT算法在不同的比特率下比EZW算法的峰值信噪比(PSNR)都有所提高，具有计算复杂度低、位速率容易控制的特点。SPIHT算法在系数子集的分割和重要信息的传输方面采用了独特的方法，能够在实现幅值大的系数优先传输的同时，隐式地传送系数的排序信息。这个隐式传送是什么意思呢？我们知道，任何排序算法的执行路径都是使用分支点的比较结果进行定义的！如果解码器和编码器使用相同的排序算法，则对于编码器输入的系数比较

5、结果，解码器通过执行相同的路径就可获得排序信息，这就是所谓的隐式传送排序信息”了。后面我们将会看到，SPIHT算法的解码、编码程序大部分代码是相同的，只在输入输出和分支点方面有所区别！二、SPIHT算法使用的树结构、分集规则和有序表1、树结构SPIHT算法的树结构与EZW算法的树结构基本相同，区别在于：对于一幅N级二维小波分解的图像，在EZW算法的零树结构中，LL_N有三个孩子HL_N、LH_N和HH_N;而SPIHT算法的树结构中，LL_N是没有孩子的！挺不好意思的说，我前面说的程序出错，就是没看清这一点，只以为是点(1,1)没有孩子，结果初始化的不重要子集表LIS就包含了具有父子关系的点，

6、造成排序扫描过程中对这些点重复扫描，生成冗余的LSP列表，重构图像失真大哎，粗心使不得啊！SPIHT算法的树结构中，树的每个节点与一个小波系数对应，我们用坐标(r,c)来标识节点或系数Cr,c。最低频子带LL_N中的系数和最高频子带中的系数没有孩子。设X是一个小波系数坐标集：X=|(r,c)|,对于正整数n,定义函数Sn(X)如下：ifmax|Cr,c|>=2AnthenSn(X)=1elseSn(X)=0如果Sn(X)=1,则坐标集X关于阈值2an是重要的，否则是不重要的。2、分集规则首先引入下面四个集合符号：(1)0(r,c)节点(r,c)所有孩子的集合；(2) D(r,c)节点(r

7、,c)所有子孙的集合(包括孩子)；(3) L(r,c)节点(r,c)所有非直系子孙的集合(即不包括孩子)；L(r,c)=D(r,c)0(r,c)(4) H所有树根的坐标集。(对N级小波分解，H就是LL_N、HL_N、LH_N和HH_N中所有系数的坐标构成的集合)根据SPIHT算法树结构的特点，除了LL_N、LL_1、HL_1、LH_1和HH_1之外，对任意系数白坐标(r,c),都有：(由于Matlab矩阵下标起始值为1,公式作了相应调整)0(r,c)=(2*r-1,2*c-1),(2*r-1,2*c),(2*r,2*c-1),(2*r,2*c)SPIHT算法的分集规则如下：(1)初始坐标集为(

8、r,c)|(r,c)CH、D(r,c)|(r,c)CH。(2)若D(r,c)关于当前阈值是重要的，则D(r,c)分裂成L(r,c)和0(r,c)。(3)若L(r,c)关于当前阈值是重要的，则L(r,c)分裂成4个集合D(r0,c0),(r0,c0)60(r,c)。3、有序表SPIHT算法引入了三个有序表来存放重要信息：(1) LSP重要系数表；(2) LIP不重要系数表；(3) LIS不重要子集表。这三个表中，每个表项都使用坐标(r,c)来标识。在LIP和LSP中，坐标(r,c)表示单个小波系数；而LIS中，坐标(r,c)代表两种系数集，即D(r,c)或L(r,c),分别称为D型表项、L型表项

9、，在Matlab实现时，用一个新的列表LisFlag来标识LIS中各个表项的类型，LisFlag的元素有D''画种字符。上一篇文章我们讨论了SPIHT算法与EZW算法的关系，介绍了SPIHT算法的树结构、分集规则和有序表的构建。在此基础上，我们接下来讨论算法的编码原理。下文给出了比较详细的数学描述，吃透了这一过程，就比较容易写出程序代码了。SPIHT算法的编码过程如下：(1)初始化输出初始阈值T的指数N=floor(log2(max|Cr,c|)(Matlab函数floor(num)给出不大于数值num的最大整数)定义：LSP为空集LIP=(r,c)|(r,c)CHLIS=D(

10、r,c)|(r,c)CH且(r,c)具有非零子孙初始的LIS中各表项类型均为口LIS和LIP中(r,c)的排列顺序与EZW算法零树结构的扫描顺序相同(即按从上到下、从左到右的“Z型次序排列)。(2)排序扫描1)扫描LIP队列对LIP队列的每个表项(r,c): 输出Sn0ut(r,c)(函数Sn0ut判断(r,c)的重要性)； 如果SnOut(r,c)=1,则向排序位流Sn输出和(r,c)的符号位(由1''0'表示)，然后将(r,c)从LIP队列中删除，添加到LSP队列的尾部。 如果SnOut(r,c)=0,则向排序位流Sn输出0'2)扫描LIS队列对LIS队列的

11、每个表项(r,c):如果(r,c)是限表项输出SnOut(D(r,c);*如果SnOut(D(r,c)=1向排序位流Sn输出1'对每个(rO,cO)6O(r,c)输出SnOut(rO,cO)*如果SnOut(rO,cO)=1,则向排序位流Sn输出和(rO,cO)的符号位，将(rO,cO)添加到LSP的尾部；*如果SnOut(rO,cO)=0,则向排序位流Sn输出0'将(rO,cO)添加到LIP的尾部。判断L(r,c)是否为空集如果L(r,c)非空，则将(r,c)作为'型表项添加到LIS的尾部；如果L(r,c)为空集，则将理表项(r,c)从LIS中删除。* 如果SnOut

12、(D(r,c)=0则向排序位流Sn输出0'如果(r,c)是'型表项输出SnOut(L(r,c);* 如果SnOut(L(r,c)=1,则向排序位流Sn输出1'然后将(r,c)的4个孩子(rO,cO)作为DS表项依次添加到LIS的尾部，将理表项(r,c)从LIS中删除；* 如果SnOut(L(r,c)=0,则向排序位流Sn输出0'(3)精细扫描将上一级扫描后的LSP列表记为LSP_Old,对于(r,c)LSP_Old,将系数Cr,c的绝对值转换为二进制表示Br,c;输出Br,c中第N个最重要的位(即对应于2AN权位处的符号或0'到精细位流Rn。(4)更新阈

13、值指数将阈值指数N减至N-1,返回到步骤(2)进行下一级编码扫描。上一篇文章已经详细介绍了SPIHT算法的编码过程，接下来有关编码和解码的部分就直接把代码写出来啦，我的代码里有详细的中文注释，基本上把程序的每个步骤都作了说明，呵呵，利人也利己！1、首先给出编码的主程序functionT,SnList,RnList,ini_LSP,ini_LIP,ini_LIS,ini_LisFlag=spihtcoding(DecIm,imDim,codeDim)%函数SPIHTCODING()是SPIHT算法的编码主程序%输入参数：DecIm小波分解系数矩阵；%imDim小波分解层数；%codeDim编码级

14、数。%输出参数：T初始阈值，T=2AN,N=floor(log2(max|c(i,j)|),c(i,j)为小波系数矩阵的元素%SnList排序扫描输出位流%RnList精细扫描输出位流%ini_L*初始系数(集合)表%LSP:重要系数表%LIP：不重要系数表%LIS:不重要子集表，其中的表项是D型或L型表项的树根点%LisFlag:LIS中各表项的类型，包括D型和L型两种globalMatrMatcMat%Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用%rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用%Threshold%Mat=DecIm

15、;MaxMat=max(max(abs(Mat);N=floor(log2(MaxMat);T=2AN;%公式：N=floor(log2(max|c(i,j)|),c(i,j)为小波系数矩阵的元素%OutputIntialization%SnList=;RnList=;ini_LSP=;ini_LIP=coef_H(imDim);rlip=size(ini_LIP,1);ini_LIS=ini_LIP(rlip/4+1:end,:);rlis=size(ini_LIS,1);ini_LisFlag(1:rlis)='D'%ini_LSP:扫描开始前无重要系数，故LSP=;%in

16、i_LIP:所有树根的坐标集，对于N层小波分解，LIP是LL_N,LH_N,HL_N,HH_N所有%系数的坐标集合；%函数COEF_H()用于计算树根坐标集H%ini_LIS:初始时，LIS是LH_N,HL_N,HH_N所有系数的坐标集合；在SPIHT算法中，%LL_N没有孩子。%ini_LisFlag:初始时，LIS列表的表项均为D型。%CodingInputInitialization%LSP=ini_LSP;LIP=ini_LIP;LIS=ini_LIS;LisFlag=ini_LisFlag;%将待输出的各项列表存入相应的编码工作列表%CodingLoop%ford=1:codeDim

17、%CodingInitialization%Sn=;LSP_Old=LSP;%每级编码产生的Sn都是独立的，故Sn初始化为空表%列表LSP_Old用于存储上级编码产生的重要系数列表LSP,作为本级精细扫描的输入%SortingPass%-LIPScan%Sn,LSP,LIP=lip_scan(Sn,N,LSP,LIP);%检查LIP表的小波系数，更新列表LIP、LSP和排序位流Sn%LISScan-%LSP,LIP,LIS,LisFlag,Sn,N=lis_scan(N,Sn,LSP,LIP,LIS,LisFlag);%这里，作为输出的N比作为输入的N少1,即out_N=in_N-1%各项输出

18、参数均作为下一编码级的输入%RefinementPass%Rn=refinement(N+1,LSP_Old);%精细扫描是在当前阈值T=2AN下，扫描上一编码级产生的LSP,故输入为(N+1,LSP_Old),%输出为精细位流Rn%OutputDataflow%SnList=SnList,Sn,7;RnList=RnList,Rn,7;%数字柞为区分符，区分不同编码级的Rn、Sn位流end编码主程序中调用到的子程序有：COEF_H():用于计算树根坐标集H,生成初始的LIP队列；LIP_SCAN():检查LIP表的各个表项是否重要，更新列表LIP、LSP和排序位流Sn;LIS_SCAN():

19、检查LIS表的各个表项是否重要，更新列表LIP、LSP、LIS、LisFlag和排序位流Sn;REFINEMENT():精细扫描编码程序，输出精细位流Rn。(1)下面是计算树根坐标集H的程序functionlp=coef_H(imDim)%函数COEF_H()根据矩阵的行列数rMat、cMat和小波分解层数imDim来计算树根坐标集H%输入参数：imDim小波分解层数，也可记作N%输出参数：lprMat*cMat矩阵经N层分解后，LL_N,LH_N,HL_N,HH_N所有系数的坐标集合globalrMatcMat%rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用

20、row=rMat/2A(imDim-1);col=cMat/2A(imDim-1);%row、col是LL_N,LH_N,HL_N,HH_N组成的矩阵的行、列数lp=listorder(row,col,1,1);%因为LIP和LIS中元素(r,c)的排列顺序与EZW零树结构的扫描顺序相同%直接调用函数LISTORDER()即可获得按EZW扫描顺序排列的LIP列表(2)这里调用了函数LISTORDER()来获取按EZW扫描顺序排列的LIP列表，以下是该函数的程序代码：functionlsorder=listorder(mr,mc,pr,pc)%函数LISTORDER()生成按理递归结构排列的坐标

21、列表%函数递归原理：对一个mr*mc的矩阵，其左上角元素的坐标为(pr,pc);首先将矩阵按田”%字型分成四个对等的子矩阵，每个子矩阵的行、列数均为mr/2、mc/2,左上角元素的坐标%从上到下、从左到右分别为(pr,pc)、(pr,pc+mc/2)、(pr+mr/2,pc)、(pr+mr/2,pc+mc/2)。%把每个子矩阵再分裂成四个矩阵，如此递归分裂下去，直至最小矩阵的行列数等于2,获取最小%矩阵的四个点的坐标值，然后逐步向上回溯，即可得到按理递归结构排列的坐标列表。lso=pr,pc;pr,pc+mc/2;pr+mr/2,pc;pr+mr/2,pc+mc/2;%列表lso是父矩阵分裂成

22、四个子矩阵后，各子矩阵左上角元素坐标的集合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

23、=lm4;ls4;endlsorder=lso;lm1;lm2;lm3;lm4;%四个子矩阵结束递归回溯到父矩阵时，列表Isorder的头四个坐标值为列表lso的元素%这四个坐标值与后面的各个子矩阵的坐标元素有重叠，故需消去%当函数输出的列表长度length(lsorder)与矩阵的元素个数mr*mc*4不相等时，%就说明有坐标重叠发生。len=length(lsorder);lsorder=lsorder(len-mr*mc*4+1:len,:);本文给出SPIHT编码的精细扫描程序，其中包括一个能够将带小数的十进制数转换为二进制表示的函数，这个转换函数可以实现任意精度的二进制转换，特别是将

24、小数部分转换为二进制表示。希望对有需要的朋友有所帮助。下一篇文章将给出SPIHT的解码程序。请关注后续文章，欢迎Email联系交流。4、精细扫描程序functionRn=refinement(N,LSP_Old)%函数REFINEMENT。为精细编码程序，对上一级编码产生白重要系数列表LSP_Old,读取每个%表项相应小波系数绝对值的二进制表示，输出其中第N个重要的位，即相应于2AN处的码数%输入参数：N本级编码阈值的指数%LSP_Old上一级编码产生的重要系数列表%输出参数：Rn精细扫描输出位流globalMat%Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用Rn=;%

25、每级精细扫描开始时，Rn均为空表%LSP_Old非空时才执行精细扫描程序ifisempty(LSP_Old)rlsp=size(LSP_Old,1);%获取LSP_Old的表项个数，对每个表项进行扫描forr=1:rlsptMat=Mat(LSP_Old(r,1),LSP_Old(r,2);%读取该表项对应的小波系数值biLSP,Np=fracnum2bin(abs(tMat),N);%函数FRACNUM2BIN()根据精细扫描对应的权位N,将任意的十进制正数转换为二进制数，%输出参数为二进制表示列表biLSP和权位N与最高权位的距离Np。Rn=Rn,biLSP(Np);%biLSP(Np)即

26、为小波系数绝对值的二进制表示中第N个重要的位endend(1)十进制数转换为二进制表示的程序functionbinlist,qLpoint=fracnum2bin(num,qLevel)%函数FRACNUM2BIN()根据精细扫描对应的权位N,将任意的十进制正数转换为二进制数，%包括带有任意位小数的十进制数。Matlab中的函数dec2bin()、dec2binvec()只能将十%进制数的整数部分转换为二进制表示，对小数部分则不转换。%输入参数：num非负的十进制数%qLevel量化转换精度，也可以是精细扫描对应的权位N%输出参数：biLSP二进制表示列表%Np权位N与最高权位的距离，N也是本

27、级编码阈值的指数intBin=dec2binvec(num);%首先用Matlab函数dec2binvec()获取整数部分的二进制表示intBin，低位在前，高位在后intBin=intBin(end:-1:1);%根据个人习惯，将二进制表示转换为高位在前，低位在后lenIB=length(intBin);%求出二进制表示的长度decpart=num-floor(num);%求出小数部分decBin=;%小数部分的二进制表示初始化为空表%根据量化精度要求输出总的二进制表示列表if(qLevel+1)>lenIB%如果量化精度高于整数部分的二进制码长，则输出为零值列表binlist=zer

28、os(1,qLevel+1);qLpoint=1;elseifqLevel>=0%如果量化精度在整数权位，则输出整数部分的二进制表示intBin%不需转换小数部分，同时输出量化精度与最高权位的距离Npbinlist=intBin;binlist(lenIB-qLevel+1:end)=0;qLpoint=lenIB-qLevel;elseifqLevel<0%如果量化精度在小数权位，则需转换小数部分N=-1;whileN>=qLevel%小数部分的转换只需进行到量化精度处res=decpart-2AN;ifres=0decBin=decBin,1;decBin(end+1:-

29、qLevel)=0;%如果小数部分的转换完成时仍未达到量化精度所在的权位，则补零break;elseifres>0decBin=decBin,1;decpart=res;N=N-1;elsedecBin=decBin,0;N=N-1;endendbinlist=intBin,decBin;qLpoint=lenIB-qLevel;%输出整数部分和小数部分的二进制表示intBin,decBin,以及量化精度与最高权位的距离Npend至此，SPIHT算法的编码程序就介绍完毕啦！以后有时间的话会增加嫡编码的功能(例如Huffman编码)。现在我们讨论SPIHT算法的解码过程。SPIHT的编码输

30、出包括：初始阈值T,排序扫描位流SnList,精细扫描位流RnList,初始有序表(LSP、LIP、LIS、LisFlag),这些参数就作为SPIHT解码的输入，另外还有解码级数decodeDim。前面我们提到，任何排序算法的执行路径都是使用分支点的比较结果进行定义的。如果解码器和编码器使用相同的排序算法，则对于编码器输入的系数比较结果，解码器通过执行相同的路径就可获得排序信息。所以，只需将编码器数学表述中的输出”改为输入"，解码器即可恢复数据的排序信息；在恢复数据排序信息的同时，解码器还要负责图像的重构，对于确认恢复的重要系数，通过排序扫描和精细扫描两个步骤更新系数的量化值，逐步提

31、高逼近精度和重构图像的质量。1、首先给出解码主程序functionDecodeMat=spihtdecoding(T,SnList,RnList,ini_LSP,ini_LIP,ini_LIS,ini_LisFlag,decodeDim)%函数SPIHTDECODING()是SPIHT算法的解码主程序%输入参数：T初始阈值，T=2AN,N=floor(log2(max|c(i,j)|),c(i,j)为小波系数矩阵的元素%SnList排序扫描输出位流%RnList精细扫描输出位流%ini_L*初始系数(集合)表%LSP：重要系数表%LIP：不重要系数表%LIS：不重要子集表，其中的表项是D型或L

32、型表项的树根点%LisFlag:LIS中各表项的类型，包括D型和L型两种%decodeDim解码级数%输出参数：DecodeMat解码后重构的小波系数矩阵%由SPIHT算法原理的排序特点，解码器的执行程序与编码器的程序路径基本相同，只需把编码器程序中的%输出信息改为输入，稍作修改即可得到解码器的程序代码globalrMatcMat%rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用%DecodingInputInitialization%N=log2(T);%获取初始阈值的指数-NDecodeMat=2A(N-decodeDim)*rand(rMat,cMat);

33、%初始化重构矩阵为一个随机矩阵，其元素最大值小于最高级解码阈值的二分之一%这样就可以保证未被扫描赋值的区域有一定的灰度，避免重构图像出现色块LSP=ini_LSP;LIP=ini_LIP;LIS=ini_LIS;LisFlag=ini_LisFlag;%将输入的各项列表存入相应的解码工作列表%DecodingLoop%ford=1:decodeDim%DecodingInitialization%Sn,SnList=getflow(SnList);Rn,RnList=getflow(RnList);%用GETFLOW()函数读取本级解码所需的位流信息LSP_Old=LSP;%列表LSP_Old

34、用于存储上级解码产生的重要系数列表LSP,作为本级精细解码的输入%SortingPass%-LIPScan-%DecodeMat,Sn,LSP,LIP=lip_decode(DecodeMat,Sn,N,LSP,LIP);%-LISScan-%LSP,LIP,LIS,LisFlag,DecodeMat,N=lis_decode(DecodeMat,N,Sn,LSP,LIP,LIS,LisFlag);%RefinementPass%DecodeMat=decRefine(DecodeMat,Rn,N+1,LSP_Old);end可以看出，解码程序和编码程序几乎是一样的，执行路径也相同。只不过解码

35、器要对输入的扫描位流SnList、RnList进行分段读取，只读入本级解码所需的位流，下面给出位流信息读取程序：(1)位流信息读取程序functionflow,bitflow=getflow(bitflow)%函数GETFLOW()用于截取本级解码所需的位流信息%输入参数：bitflow初始为编码器的输入位流，在解码过程中为上一级解码截取后剩余的编码位流%输出参数：flow本级解码所需的位流(排序位流Sn、精细位流Rn)%bitflow本级解码截取后剩余的编码位流flow=;i=1;whilebitflow(i)=7flow(i)=bitflow(i);i=i+1;end%数字作为位流区别符，

36、区分不同编码级的位流bitflow(1:i)=;%将对应于本级解码位流的bitflow数据清空%bitflow被完全清空后返回一个空表值ifisempty(bitflow)bitflow=;end提示：任何排序算法的执行路径都是使用分支点的比较结果进行定义的。如果解码器和编码器使用相同的排序算法，则对于编码器输入的系数比较结果，解码器通过执行相同的路径就可获得排序信息。所以，只需将编码器数学表述中的输出”改为输入"，解码器即可恢复数据的排序信息；在恢复数据排序信息的同时，解码器还要负责图像的重构，对于确认恢复的重要系数，通过排序扫描和精细扫描两个步骤更新系数的量化值，逐步提高逼近精度和重构图像的质量。globalrMatcMat%rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相