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

忘是亡心i 2021-12-12 01:17 330阅读 0赞

本文来自CSDN博客，：[http://blog.csdn.net/chenyusiyuan/archive/2007/12/07/1923334.aspx][http_blog.csdn.net_chenyusiyuan_archive_2007_12_07_1923334.aspx]

多级树集合分裂（SPIHT）算法的过程详解与Matlab实现（1）算法概述 收藏  
   
上星期我们讨论了EZW算法，很高兴收到了一些朋友的email，对算法进行探讨、交流。这也是我开这个博客的源动力之一，学习就应该开诚布公、交流互助，在探讨中加深对所学知识的理解和掌握。在弄懂了EZW算法原理并用Matlab实现后，我继续学习EZW的改进算法。至今有一周的时间没更新博客、写新文章了，其实就是把时间用在EZW的一个改进算法——多级树集合分裂（Set Partitioning in Hierarchical Trees, SPIHT）算法的研究和Matlab实现。由于EZW是SPIHT的基础，所以在EZW算法的Matlab代码的基础上，我很快就完成了SPIHT的代码编写，但最痛苦的是一开始没吃透算法原理，程序在初始分集上出了错，调试了两天找不出根本问题，昨天从头再看一次算法原理，才发现问题所在……呵呵，小小的粗心就耽搁了我两三天的时间和精力！问题解决后，就编写程序注释了，上次EZW算法的代码都没写注释，让大家看着辛苦，不好意思哦！好，接下来就开始讨论SPIHT算法的原理，然后给出具体的Matlab代码。

一、SPIHT算法与EZW算法  
EZW算法是一种基于零树的嵌入式图象编码算法，虽然在小波变换系数中，零树是一个比较有效的表示不重要系数的数据结构，但是，在小波系数中还存在这样的树结构，它的树根是重要的，除树根以外的其它结点是不重要的。对这样的系数结构，零树就不是一种很有效的表示方法。A.Said和W.A.Pearlman根据Shapiro零树编码算法(EZW)的基本思想，提出了一种新的且性能更优的实现方法，即基于多级树集合分裂排序（Set Partitioning in Hierarchical Trees, SPIHT）的编码算法。它采用了空间方向树（SOT：spatial orientation tree）、全体子孙集合D（i,j）和非直系子孙集合L(i,j)的概念以更有效地表示上述特征的系数结构，从而提高了编码效率。  
SPIHT算法能够生成一个嵌入位流（embedded bit stream），使接收的位流在任意点中断时，都可解压和重构图像，具有良好的渐进传输特性；算法的初始化过程、细化过程类似于EZW算法，它改进了EZW 重要图的表示方法，也就是重要系数在表中的排序信息，使得集合的表示更为精简，从而提高了编码效率和图像压缩率。SPIHT算法在不同的比特率下比EZW算法的峰值信噪比（PSNR）都有所提高，具有计算复杂度低、位速率容易控制的特点。

SPIHT算法在系数子集的分割和重要信息的传输方面采用了独特的方法，能够在实现幅值大的系数优先传输的同时，隐式地传送系数的排序信息。这个隐式传送是什么意思呢？我们知道，任何排序算法的执行路径都是使用分支点的比较结果进行定义的！如果解码器和编码器使用相同的排序算法，则对于编码器输入的系数比较结果，解码器通过执行相同的路径就可获得排序信息，这就是所谓的“隐式传送排序信息”了。后面我们将会看到，SPIHT算法的解码、编码程序大部分代码是相同的，只在输入输出和分支点方面有所区别！  
   
二、SPIHT算法使用的树结构、分集规则和有序表  
1、树结构  
SPIHT算法的树结构与EZW算法的树结构基本相同，区别在于：  
对于一幅N级二维小波分解的图像，在EZW算法的零树结构中，LL\_N有三个孩子HL\_N、LH\_N和HH\_N；而SPIHT算法的树结构中，LL\_N是没有孩子的！  
挺不好意思的说，我前面说的程序出错，就是没看清这一点，只以为是点（1，1）没有孩子，结果初始化的不重要子集表LIS就包含了具有父子关系的点，造成排序扫描过程中对这些点重复扫描，生成冗余的LSP列表，重构图像失真大……哎，粗心使不得啊！  
SPIHT算法的树结构中，树的每个节点与一个小波系数对应，我们用坐标（r,c）来标识节点或系数Cr,c。最低频子带LL\_N中的系数和最高频子带中的系数没有孩子。  
设X是一个小波系数坐标集：X=\{| (r,c) |\}，对于正整数n，定义函数Sn (X) 如下：  
if max\{| Cr,c |\}>= 2 ^ n then Sn (X) = 1

else Sn (X) = 0

如果Sn (X) = 1，则坐标集X关于阈值2 ^ n 是重要的，否则是不重要的。  
   
2、分集规则  
首先引入下面四个集合符号：  
      （1）O (r,c) —— 节点(r,c)所有孩子的集合；  
      （2）D (r,c) —— 节点(r,c)所有子孙的集合（包括孩子）；  
      （3）L (r,c) —— 节点(r,c)所有非直系子孙的集合（即不包括孩子）；  
L (r,c) = D (r,c) — O (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，公式作了相应调整）  
O (r,c) = \{ (2\*r-1,2\*c-1), (2\*r-1,2\*c), (2\*r,2\*c-1), (2\*r,2\*c) \}

SPIHT算法的分集规则如下：  
      （1） 初始坐标集为\{(r,c) | (r,c)∈H \}、\{D(r,c) | (r,c)∈H \}。  
      （2） 若D(r,c) 关于当前阈值是重要的，则D(r,c) 分裂成 L (r,c) 和O (r,c)。  
      （3） 若L (r,c) 关于当前阈值是重要的，则L (r,c) 分裂成 4 个集合 D (rO,cO)， (rO,cO) ∈O (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型表项，在Matlab实现时，用一个新的列表 LisFlag 来标识 LIS 中各个表项的类型，LisFlag的元素有‘D’、‘L’两种字符。

上一篇文章我们讨论了SPIHT算法与EZW算法的关系，介绍了SPIHT算法的树结构、分集规则和有序表的构建。在此基础上，我们接下来讨论算法的编码原理。下文给出了比较详细的数学描述，吃透了这一过程，就比较容易写出程序代码了。

SPIHT算法的编码过程如下：  
（1）初始化  
      输出初始阈值T的指数 N = floor ( log2 ( max\{| Cr,c |\} ) ) （Matlab函数 floor( num ) 给出不大于数值 num 的最大整数）  
      定义：  LSP 为空集  
                  LIP = \{(r,c) | (r,c)∈H \}  
                  LIS = \{D(r,c) | (r,c)∈H 且(r,c)具有非零子孙\}  
      初始的LIS中各表项类型均为‘D’， LIS 和 LIP 中 (r,c) 的排列顺序与EZW算法零树结构的扫描顺序相同（即按从上到下、从左到右的“Z”型次序排列）。  
（2）排序扫描  
      1）扫描LIP队列  
      对LIP队列的每个表项 (r,c) ：  
            ①    输出SnOut（r,c）（函数SnOut判断(r,c)的重要性）；  
            ②    如果SnOut（r,c）= 1，则向排序位流Sn输出‘1’和（r,c）的符号位（由‘1’、‘0’表示），然后将（r,c）从LIP队列中删除，添加到LSP队列的尾部。  
            ③    如果SnOut（r,c）= 0，则向排序位流Sn输出‘0’。  
      2）扫描LIS队列  
      对LIS队列的每个表项 (r,c) ：  
            ①    如果（r,c）是‘D’型表项  
                  输出SnOut（D (r,c)）；  
                  \* 如果SnOut（D (r,c)）= 1  
                          向排序位流Sn输出‘1’；  
                          对每个(rO,cO) ∈O (r,c)  
                              输出SnOut（rO,cO）  
                              \* 如果SnOut（rO,cO）= 1，则向排序位流Sn输出‘1’和（rO,cO）的符号位，将（rO,cO）添加到LSP的尾部；  
                              \* 如果SnOut（rO,cO）= 0，则向排序位流Sn输出‘0’，将（rO,cO）添加到LIP的尾部。  
                          判断L (r,c) 是否为空集  
                              如果L (r,c) 非空，则将（r,c）作为‘L’型表项添加到LIS的尾部；  
                              如果L (r,c)为空集，则将‘D’型表项（r,c）从LIS中删除。  
                  \* 如果SnOut（D (r,c)）= 0  
                          则向排序位流Sn输出‘0’。  
            ②    如果（r,c）是‘L’型表项  
                  输出SnOut（L (r,c)）；  
                  \* 如果SnOut（L (r,c)）= 1，则向排序位流Sn输出‘1’，然后将(r,c)的4个孩子（rO,cO）作为‘D’型表项依次添加到LIS的尾部，将‘L’型表项（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个最重要的位（即对应于2^N权位处的符号‘1’或‘0’）到精细位流Rn。  
（4）更新阈值指数  
      将阈值指数N减至N—1，返回到步骤（2）进行下一级编码扫描。

[ ][Link 1]

上一篇文章已经详细介绍了SPIHT算法的编码过程，接下来有关编码和解码的部分就直接把代码写出来啦，我的代码里有详细的中文注释，基本上把程序的每个步骤都作了说明，呵呵，利人也利己！

1、首先给出编码的主程序

function \[T,SnList,RnList,ini\_LSP,ini\_LIP,ini\_LIS,ini\_LisFlag\]=spihtcoding(DecIm,imDim,codeDim)  
% 函数 SPIHTCODING() 是SPIHT算法的编码主程序  
% 输入参数：DecIm   ——小波分解系数矩阵；  
%                     imDim   ——小波分解层数；  
%                     codeDim ——编码级数。  
% 输出参数：T —— 初始阈值，T=2^N，N=floor(log2(max\{|c(i,j)|\}))，c(i,j)为小波系数矩阵的元素  
%                     SnList —— 排序扫描输出位流  
%                     RnList —— 精细扫描输出位流  
%                     ini\_L\* —— 初始系数（集合）表  
%                     LSP：重要系数表  
%                     LIP：不重要系数表  
%                     LIS：不重要子集表，其中的表项是D型或L型表项的树根点  
%                     LisFlag：LIS中各表项的类型，包括D型和L型两种

global Mat rMat cMat  
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用  
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用

%---------------------------%  
% ----- Threshold ----- %  
%---------------------------%  
Mat=DecIm;  
MaxMat=max(max(abs(Mat)));  
N=floor(log2(MaxMat));  
T=2^N;  
% 公式：N=floor(log2(max\{|c(i,j)|\}))，c(i,j)为小波系数矩阵的元素

%--------------------------------------%  
% ----- Output Intialization ----- %  
%--------------------------------------%  
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=\[\]；  
% ini\_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型。

%------------------------------------------------%  
% ----- Coding Input Initialization ------ %  
%------------------------------------------------%  
LSP=ini\_LSP;  
LIP=ini\_LIP;  
LIS=ini\_LIS;  
LisFlag=ini\_LisFlag;  
% 将待输出的各项列表存入相应的编码工作列表

%--------------------------------%  
% ----- Coding Loop ------ %  
%--------------------------------%  
for d=1:codeDim  
    %------------------------------------------%  
    % ----- Coding Initialization ------- %  
    %------------------------------------------%  
    Sn=\[\];  
    LSP\_Old=LSP;  
    % 每级编码产生的Sn都是独立的，故Sn初始化为空表  
    % 列表LSP\_Old用于存储上级编码产生的重要系数列表LSP，作为本级精细扫描的输入  
    %-------------------------------%  
    % ----- Sorting Pass ----- %  
    %-------------------------------%  
        % ----- LIP Scan -------- %  
        %----------------------------%  
    \[Sn,LSP,LIP\]=lip\_scan(Sn,N,LSP,LIP);  
    % 检查LIP表的小波系数，更新列表LIP、LSP和排序位流 Sn  
        %-------------------------%  
        % ----- LIS Scan ----- %  
        %-------------------------%  
    \[LSP,LIP,LIS,LisFlag,Sn,N\]=lis\_scan(N,Sn,LSP,LIP,LIS,LisFlag);  
    % 这里，作为输出的N比作为输入的N少1，即 out\_N=in\_N-1  
    % 各项输出参数均作为下一编码级的输入  
    %-------------------------------------%  
    % ----- Refinement Pass ----- %  
    %-------------------------------------%  
    Rn=refinement(N+1,LSP\_Old);    
    % 精细扫描是在当前阈值T=2^N下，扫描上一编码级产生的LSP，故输入为(N+1,LSP\_Old)，  
    % 输出为精细位流 Rn  
    %-----------------------------------%  
    % ----- Output Dataflow ----- %  
    %-----------------------------------%  
    SnList=\[SnList,Sn,7\];  
    RnList=\[RnList,Rn,7\];  
    % 数字‘7’作为区分符,区分不同编码级的Rn、Sn位流  
end

编码主程序中调用到的子程序有：

COEF\_H() ：用于计算树根坐标集 H ，生成初始的LIP队列；  
LIP\_SCAN() ：检查LIP表的各个表项是否重要，更新列表LIP、LSP和排序位流 Sn ；  
LIS\_SCAN() ：检查LIS表的各个表项是否重要，更新列表LIP、LSP、LIS、LisFlag和排序位流 Sn ；  
REFINEMENT() ：精细扫描编码程序，输出精细位流 Rn 。

（1）下面是计算树根坐标集 H 的程序

function lp=coef\_H(imDim)  
% 函数 COEF\_H() 根据矩阵的行列数rMat、cMat和小波分解层数imDim来计算树根坐标集 H  
% 输入参数：imDim —— 小波分解层数,也可记作 N  
% 输出参数：lp —— rMat\*cMat矩阵经N层分解后，LL\_N,LH\_N,HL\_N,HH\_N 所有系数的坐标集合

global rMat cMat  
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用

row=rMat/2^(imDim-1);  
col=cMat/2^(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列表，以下是该函数的程序代码：

function lsorder=listorder(mr,mc,pr,pc)  
% 函数 LISTORDER() 生成按‘Z’型递归结构排列的坐标列表  
% 函数递归原理：对一个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，获取最小  
% 矩阵的四个点的坐标值，然后逐步向上回溯，即可得到按‘Z’型递归结构排列的坐标列表。

lso=\[pr,pc;pr,pc+mc/2;pr+mr/2,pc;pr+mr/2,pc+mc/2\];  
% 列表lso是父矩阵分裂成四个子矩阵后，各子矩阵左上角元素坐标的集合  
mr=mr/2;  
mc=mc/2;  
% 子矩阵的行列数是父矩阵的一半  
lm1=\[\];lm2=\[\];lm3=\[\];lm4=\[\];  
if (mr>1)&&(mc>1)  
% 按‘Z’型结构递归  
    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\];  
% 四个子矩阵结束递归回溯到父矩阵时，列表lsorder的头四个坐标值为列表lso的元素  
% 这四个坐标值与后面的各个子矩阵的坐标元素有重叠，故需消去  
% 当函数输出的列表长度length(lsorder)与矩阵的元素个数mr\*mc\*4不相等时，  
% 就说明有坐标重叠发生。  
len=length(lsorder);  
lsorder=lsorder(len-mr\*mc\*4+1:len,:);

本文给出SPIHT编码的排序扫描代码，排序扫描分为LIP队列扫描和LIS队列扫描两个步骤，其中LIS队列扫描较为复杂，在编程时容易出现错误，要倍加注意。

2、LIP队列扫描程序

function \[Sn,LSP,LIP\]=lip\_scan(Sn,N,LSP,LIP)  
% 函数 LIP\_SCAN() 检查LIP表的各个表项是否重要，更新列表LIP、LSP和排序位流 Sn  
% 输入参数：Sn —— 本级编码排序位流，为空表  
%                     N —— 本级编码阈值的指数  
%                     LSP —— 上一级编码生成的重要系数列表  
%                     LIP —— 上一级编码生成的不重要系数列表  
% 输出参数：Sn —— 对上一级编码生成的LIP列表扫描后更新的排序位流  
%                     LSP —— 对上一级编码生成的LIP列表扫描后更新的重要系数列表  
%                     LIP —— 经本级LIP扫描处理后更新的不重要系数列表

global Mat  
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用

rlip=size(LIP,1);  
% r 是指向 LIP 当前读入表项位置的指针  
r=1;  
% 由于循环过程中列表 LIP 的表长会变化，不适合用 for 循环，故采用 while 循环  
while r<=rlip  
% 读入当前表项的坐标值  
    rN=LIP(r,1);  
    cN=LIP(r,2);  
% 调用 SNOUT() 函数来判断该表项是否重要  
    if SnOut(LIP(r,:),N)  
        % 若重要，则输入‘1’到 Sn  
        Sn=\[Sn,1\];  
        % 输入正负符号‘1’或‘0’到 Sn  
        if Mat(rN,cN)>=0  
            Sn=\[Sn,1\];  
        else  
            Sn=\[Sn,0\];  
        end  
        % 将该表项添加到重要系数列表 LSP  
        LSP=\[LSP;LIP(r,:)\];  
        % 将该表项从 LIP 中删除  
        LIP(r,:)=\[\];  
    else  
        % 若不重要，则输入‘0’到 Sn  
        Sn=\[Sn,0\];  
        % 将指针指向下一个表项  
        r=r+1;  
    end  
    % 判断当前 LIP 的表长  
    rlip=size(LIP,1);  
end

3、LIS队列扫描程序

function \[LSP,LIP,LIS,LisFlag,Sn,N\]=lis\_scan(N,Sn,LSP,LIP,LIS,LisFlag)  
% 函数 LIS\_SCAN() 检查LIS表的各个表项是否重要，更新列表LIP、LSP、LIS、LisFlag和排序位流 Sn  
% 输入参数：N —— 本级编码阈值的指数  
%                     Sn —— 本级编码排序位流，为空表  
%                     LSP —— 上一级编码生成的重要系数列表  
%                     LIP —— 上一级编码生成的不重要系数列表  
%                     LIS —— 上一级编码生成的不重要子集列表  
%                     LisFlag —— 上一级编码生成的不重要子集表项类型列表  
% 输出参数：LSP —— 本级LIS列表扫描后更新的重要系数列表  
%                     LIP —— 经本级LIS扫描处理后更新的不重要系数列表  
%                     LIS —— 本级LIS列表扫描后更新的不重要子集列表  
%                     LisFlag —— 本级LIS列表扫描后更新的不重要子集表项类型列表  
%                     Sn —— 本级LIS列表扫描后更新的排序位流  
%                     N —— 下一级编码阈值的指数

% 读入当前 LIS 的表长  
rlis=size(LIS,1);  
% ls 是指向 LIS 当前表项位置的指针，初始位置为1  
ls=1;  
% 由于循环过程中列表 LIS 的表长会变化，不适合用 for 循环，故采用 while 循环  
while ls<=rlis  
    % 读入当前 LIS 表项的类型  
    switch LisFlag(ls)  
        % ‘D’类表项，包含孩子和非直系子孙  
        case 'D'  
            % 读入该表项的坐标值  
            rP=LIS(ls,1);  
            cP=LIS(ls,2);  
            % 生成‘D’型子孙树  
            chD=coef\_DOL(rP,cP,'D');  
            % 判断该子孙树是否重要  
            isImt=SnOut(chD,N);  
            if isImt  
                % 如果该子孙树重要，就输入‘1’到 Sn  
                Sn=\[Sn,1\];  
                % 生成该表项的孩子树  
                chO=coef\_DOL(rP,cP,'O');  
                % 分别判断每个孩子的重要性  
                for r=1:4  
                    % 判断该孩子的重要性和正负符号  
                    \[isImt,Sign\]=SnOut(chO(r,:),N);  
                    if isImt  
                        % 如果重要，就输入‘1’和正负标志到 Sn  
                        Sn=\[Sn,1\];  
                        if Sign  
                            Sn=\[Sn,1\];  
                        else  
                            Sn=\[Sn,0\];  
                        end  
                        % 将该孩子添加到重要系数列表 LSP  
                        LSP=\[LSP;chO(r,:)\];  
                    else  
                        % 如果不重要，就输入‘0’到 Sn  
                        Sn=\[Sn,0\];  
                        % 将该孩子添加到不重要系数列表 LIP  
                        % 本级阈值下不重要的系数在下一级编码中可能是重要的  
                        LIP=\[LIP;chO(r,:)\];  
                    end  
                end  
                % 生成该表项的非直系子孙树  
                chL=coef\_DOL(rP,cP,'L');  
                if ~isempty(chL)  
                    % 如果‘L’型树非空，则将该表项添加到列表 LIS 的尾端等待扫描  
                    LIS=\[LIS;LIS(ls,:)\];  
                    % 表项类型更改为‘L’型  
                    LisFlag=\[LisFlag,'L'\];  
                    % 至此，该表项的‘D’型LIS扫描结束，在LIS中删除该项及其类型符  
                    LIS(ls,:)=\[\];  
                    LisFlag(ls)=\[\];  
                else  
                    % 如果‘L’型树为空集  
                    % 则该表项的‘D’型LIS扫描结束，在LIS中删除该项及其类型符  
                    LIS(ls,:)=\[\];  
                    LisFlag(ls)=\[\];  
                end  
            else  
                % 如果该表项的‘D’型子孙树不重要，就输入‘0’到 Sn  
                Sn=\[Sn,0\];  
                % 将指针指向下一个 LIS 表项  
                ls=ls+1;  
            end  
            % 更新当前 LIS 的表长，转入下一表项的扫描  
            rlis=size(LIS,1);  
        case 'L'  
            % 对‘L’类表项，不需判断孩子的重要性  
            % 读入该表项的坐标值  
            rP=LIS(ls,1);  
            cP=LIS(ls,2);  
            % 生成‘L’型子孙树  
            chL=coef\_DOL(rP,cP,'L');  
            % 判断该子孙树是否重要  
            isImt=SnOut(chL,N);  
            if isImt  
                % 如果该子孙树重要，就输入‘1’到 Sn  
                Sn=\[Sn,1\];  
                % 生成该表项的孩子树  
                chO=coef\_DOL(rP,cP,'O');  
                % 将该‘L’类表项从 LIS 中删除  
                LIS(ls,:)=\[\];  
                LisFlag(ls)=\[\];  
                % 将表项的四个孩子添加到 LIS 的结尾，标记为‘D’类表项  
                LIS=\[LIS;chO(1:4,:)\];  
                LisFlag(end+1:end+4)='D';  
            else  
                % 如果该子孙树是不重要的，就输入‘0’到 Sn  
                Sn=\[Sn,0\];  
                % 将指针指向下一个 LIS 表项  
                ls=ls+1;  
            end  
            % 更新当前 LIS 的表长，转入下一表项的扫描  
            rlis=size(LIS,1);  
    end  
end  
% 对 LIS 的扫描结束，将本级阈值的指数减1，准备进入下一级编码  
N=N-1;

LIS队列扫描程序中调用的子程序有：

COEF\_DOL() ：根据子孙树的类型'type'来计算点(r,c)的子孙列表；  
SNOUT() ：根据本级阈值指数 N 判断坐标集 coefSet 是否重要；

（1）子孙树生成程序

function chList=coef\_DOL(r,c,type)  
% 函数 COEF\_DOL() 根据子孙树的类型'type'来计算点(r,c)的子孙列表  
% 输入参数：r,c —— 小波系数的坐标值  
%                     type —— 子孙树的类型  
%                     'D'：点(r,c)的所有子孙（包括孩子）  
%                     'O'：点(r,c)的所有孩子  
%                     'L'：点(r,c)的所有非直系子孙，L(r,c)=D(r,c)-O(r,c)  
% 输出参数：chList —— 点(r,c)的'type'型子孙列表

\[Dch,Och\]=childMat(r,c);  
% 函数 CHILDMAT() 返回点(r,c)的'D'型和'O'型子孙列表  
Lch=setdiff(Dch,Och,'rows');  
% 根据关系式 L(r,c)=D(r,c)-O(r,c)求出'L'型子孙列表  
% Matlab函数 SETDIFF(A,B)计算具有相同列数的两个矩阵A、B中，A有而B无的元素集合  
% 根据输入参数'type'选择输出参数  
switch type  
    case 'D'  
        chList=Dch;  
    case 'O'  
        chList=Och;  
    case 'L'  
        chList=Lch;  
end

function \[trAll,trChl\]=childMat(trRows,trCols)  
% 函数 CHILDMAT() 根据输入的坐标值trRows、trCols 输出其全体子孙 trAll，  
% 其中包括孩子树 trChl；另外，根据算法原理，还要判断子孙树是否全为零，  
% 若为全零，则trAll、trChl均为空表

trAll=treeMat(trRows,trCols);  
% 调用函数 treeMat() 生成该点的子孙树坐标队列  
trZero=1;  
% 用变量 trZero 来标记该点是否具有非零子孙  
rA=size(trAll,1);  
% 如果子孙树 trAll 中有系数值不为零，则 trZero=0，表示该点具有非零子孙  
for r=1:rA  
    if Mat(trAll(r,1),trAll(r,2))~=0  
        trZero=0;  
        break;  
    end  
end  
if trZero  
    trAll=\[\];  
    trChl=\[\];  
else  
    trChl=trAll(1:4,:);  
    % 函数 treeMat() 输出的全体子孙树 trAll 头四位元素就组成相应的孩子树  
end

上面调用的函数treeMat() 与EZW算法中使用的程序是一样的，这里就不写出来了，详细代码请参见《嵌入式小波零树（EZW）算法的过程详解和Matlab代码（1）构建扫描次序表》。

（2）系数集重要性判别程序

function \[isImt,Sign\]=SnOut(coefSet,N)  
% 函数 SNOUT() 根据本级阈值指数 N 判断坐标集 coefSet 是否重要 isImt ，对单元素  
% 的系数集输出该元素的正负符号 Sign 。

global Mat  
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用  
allMat=\[\];  
isImt=0;  
Sign=0;  
% 默认坐标集是不重要的，且首位元素是负值  
rSet=size(coefSet,1);  
% 读取坐标集中各元素的系数值  
for r=1:rSet  
    allMat(r)=Mat(coefSet(r,1),coefSet(r,2));  
    if abs(allMat(r))>=2^N  
        isImt=1;  
        break;  
    end  
end  
% 对单个元素的坐标集，判断该元素的正负符号  
% 由于函数 childMat() 对子孙全零的点会返回空表，所以要检查allMat是否为空  
if ~isempty(allMat)&&(allMat(1)>=0)  
    Sign=1;  
end

本文给出SPIHT编码的精细扫描程序，其中包括一个能够将带小数的十进制数转换为二进制表示的函数，这个转换函数可以实现任意精度的二进制转换，特别是将小数部分转换为二进制表示。希望对有需要的朋友有所帮助。下一篇文章将给出SPIHT的解码程序。请关注后续文章，欢迎 Email 联系交流。

4、精细扫描程序

function Rn=refinement(N,LSP\_Old)  
% 函数 REFINEMENT()为精细编码程序，对上一级编码产生的重要系数列表LSP\_Old，读取每个  
% 表项相应小波系数绝对值的二进制表示，输出其中第N个重要的位，即相应于 2^N 处的码数  
% 输入参数：N —— 本级编码阈值的指数  
%                 LSP\_Old —— 上一级编码产生的重要系数列表  
% 输出参数：Rn —— 精细扫描输出位流

global Mat  
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用

Rn=\[\];  
% 每级精细扫描开始时，Rn 均为空表  
% LSP\_Old 非空时才执行精细扫描程序  
if ~isempty(LSP\_Old)  
    rlsp=size(LSP\_Old,1);  
% 获取 LSP\_Old 的表项个数，对每个表项进行扫描  
    for r=1:rlsp  
        tMat=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)即为小波系数绝对值的二进制表示中第N个重要的位  
    end  
end

（1）十进制数转换为二进制表示的程序

function \[binlist,qLpoint\]=fracnum2bin(num,qLevel)  
% 函数 FRACNUM2BIN() 根据精细扫描对应的权位 N ，将任意的十进制正数转换为二进制数，  
% 包括带有任意位小数的十进制数。Matlab中的函数 dec2bin()、dec2binvec()只能将十  
% 进制数的整数部分转换为二进制表示，对小数部分则不转换。  
%  
% 输入参数：num —— 非负的十进制数  
%                 qLevel —— 量化转换精度，也可以是精细扫描对应的权位 N  
% 输出参数：biLSP —— 二进制表示列表  
%                 Np —— 权位N与最高权位的距离，N 也是本级编码阈值的指数

intBin=dec2binvec(num);  
% 首先用Matlab函数dec2binvec()获取整数部分的二进制表示intBin，低位在前，高位在后  
intBin=intBin(end:-1:1);  
% 根据个人习惯，将二进制表示转换为高位在前，低位在后  
lenIB=length(intBin);  
% 求出二进制表示的长度  
decpart=num-floor(num);  
% 求出小数部分  
decBin=\[\];  
% 小数部分的二进制表示初始化为空表

% 根据量化精度要求输出总的二进制表示列表  
if (qLevel+1)>lenIB  
% 如果量化精度高于整数部分的二进制码长，则输出为零值列表  
    binlist=zeros(1,qLevel+1);  
    qLpoint=1;  
elseif qLevel>=0  
% 如果量化精度在整数权位，则输出整数部分的二进制表示intBin  
% 不需转换小数部分，同时输出量化精度与最高权位的距离Np  
    binlist=intBin;  
    binlist(lenIB-qLevel+1:end)=0;  
    qLpoint=lenIB-qLevel;  
elseif qLevel<0  
% 如果量化精度在小数权位，则需转换小数部分  
    N=-1;  
    while N>=qLevel  
    % 小数部分的转换只需进行到量化精度处  
        res=decpart-2^N;  
        if res==0  
            decBin=\[decBin,1\];  
            decBin(end+1:-qLevel)=0;  
            % 如果小数部分的转换完成时仍未达到量化精度所在的权位，则补零  
            break;  
        elseif res>0  
            decBin=\[decBin,1\];  
            decpart=res;  
            N=N-1;  
        else  
            decBin=\[decBin,0\];  
            N=N-1;  
        end  
    end  
    binlist=\[intBin,decBin\];  
    qLpoint=lenIB-qLevel;  
    % 输出整数部分和小数部分的二进制表示intBin,decBin，以及量化精度与最高权位的距离Np  
end  
     
至此，SPIHT算法的编码程序就介绍完毕啦！以后有时间的话会增加熵编码的功能（例如Huffman编码）。

现在我们讨论SPIHT算法的解码过程。SPIHT的编码输出包括：初始阈值T，排序扫描位流SnList，精细扫描位流 RnList，初始有序表（LSP、LIP、LIS、LisFlag），这些参数就作为SPIHT解码的输入，另外还有解码级数decodeDim。

前面我们提到，任何排序算法的执行路径都是使用分支点的比较结果进行定义的。如果解码器和编码器使用相同的排序算法，则对于编码器输入的系数比较结果，解码器通过执行相同的路径就可获得排序信息。所以，只需将编码器数学表述中的“输出”改为“输入”，解码器即可恢复数据的排序信息；在恢复数据排序信息的同时，解码器还要负责图像的重构，对于确认恢复的重要系数，通过排序扫描和精细扫描两个步骤更新系数的量化值，逐步提高逼近精度和重构图像的质量。

1、首先给出解码主程序

function DecodeMat=spihtdecoding(T,SnList,RnList,ini\_LSP,ini\_LIP,ini\_LIS,ini\_LisFlag,decodeDim)  
% 函数 SPIHTDECODING() 是SPIHT算法的解码主程序  
% 输入参数：T —— 初始阈值，T=2^N，N=floor(log2(max\{|c(i,j)|\}))，c(i,j)为小波系数矩阵的元素  
%          SnList —— 排序扫描输出位流  
%          RnList —— 精细扫描输出位流  
%          ini\_L\* —— 初始系数（集合）表  
%               LSP：重要系数表  
%               LIP：不重要系数表  
%               LIS：不重要子集表，其中的表项是D型或L型表项的树根点  
%               LisFlag：LIS中各表项的类型，包括D型和L型两种  
%          decodeDim —— 解码级数  
% 输出参数：DecodeMat —— 解码后重构的小波系数矩阵  
%  
% 由SPIHT算法原理的排序特点，解码器的执行程序与编码器的程序路径基本相同，只需把编码器程序中的  
% 输出信息改为输入，稍作修改即可得到解码器的程序代码

global rMat cMat  
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用

%-------------------------------------------%  
% ----- Decoding Input Initialization ----- %  
%-------------------------------------------%  
N=log2(T);  
% 获取初始阈值的指数-N  
DecodeMat=2^(N-decodeDim)\*rand(rMat,cMat);  
% 初始化重构矩阵为一个随机矩阵，其元素最大值小于最高级解码阈值的二分之一  
% 这样就可以保证未被扫描赋值的区域有一定的灰度，避免重构图像出现色块  
LSP=ini\_LSP;  
LIP=ini\_LIP;  
LIS=ini\_LIS;  
LisFlag=ini\_LisFlag;  
% 将输入的各项列表存入相应的解码工作列表

%---------------------------%  
% ----- Decoding Loop ----- %  
%---------------------------%  
for d=1:decodeDim  
    %-------------------------------------%  
    % ----- Decoding Initialization ----- %  
    %-------------------------------------%  
    \[Sn,SnList\]=getflow(SnList);  
    \[Rn,RnList\]=getflow(RnList);  
% 用 GETFLOW() 函数读取本级解码所需的位流信息  
    LSP\_Old=LSP;  
% 列表LSP\_Old用于存储上级解码产生的重要系数列表LSP，作为本级精细解码的输入

%--------------------------%  
    % ----- Sorting Pass ----- %  
    %--------------------------%  
        % ----- LIP Scan ----- %  
        %----------------------%  
    \[DecodeMat,Sn,LSP,LIP\]=lip\_decode(DecodeMat,Sn,N,LSP,LIP);  
        %----------------------%  
        % ----- LIS Scan ----- %  
        %----------------------%  
    \[LSP,LIP,LIS,LisFlag,DecodeMat,N\]=lis\_decode(DecodeMat,N,Sn,LSP,LIP,LIS,LisFlag);  
    %-----------------------------%  
    % ----- Refinement Pass ----- %  
    %-----------------------------%  
    DecodeMat=decRefine(DecodeMat,Rn,N+1,LSP\_Old);  
end

可以看出，解码程序和编码程序几乎是一样的，执行路径也相同。只不过解码器要对输入的扫描位流SnList、RnList进行分段读取，只读入本级解码所需的位流，下面给出位流信息读取程序：

（1）位流信息读取程序

function \[flow,bitflow\]=getflow(bitflow)  
% 函数 GETFLOW() 用于截取本级解码所需的位流信息  
% 输入参数：bitflow —— 初始为编码器的输入位流，在解码过程中为上一级解码截取后剩余的编码位流  
% 输出参数：flow —— 本级解码所需的位流（排序位流Sn、精细位流Rn）  
%                 bitflow —— 本级解码截取后剩余的编码位流

flow=\[\];  
i=1;  
while bitflow(i)~=7  
    flow(i)=bitflow(i);  
    i=i+1;  
end  
% 数字‘7’作为位流区别符，区分不同编码级的位流  
bitflow(1:i)=\[\];  
% 将对应于本级解码位流的 bitflow 数据清空  
% bitflow 被完全清空后返回一个空表值  
if isempty(bitflow)  
    bitflow=\[\];  
end

提示：任何排序算法的执行路径都是使用分支点的比较结果进行定义的。如果解码器和编码器使用相同的排序算法，则对于编码器输入的系数比较结果，解码器通过执行相同的路径就可获得排序信息。

所以，只需将编码器数学表述中的“输出”改为“输入”，解码器即可恢复数据的排序信息；在恢复数据排序信息的同时，解码器还要负责图像的重构，对于确认恢复的重要系数，通过排序扫描和精细扫描两个步骤更新系数的量化值，逐步提高逼近精度和重构图像的质量。

global rMat cMat  
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用

% 读入当前 LIS 的表长  
rlis=size(LIS,1);  
% ls 是指向 LIS 当前表项位置的指针，初始位置为1  
ls=1;  
while ls<=rlis  
    % 读入当前 LIS 表项的类型  
    switch LisFlag(ls)  
        % ‘D’类表项，包含孩子和非直系子孙  
        case 'D'  
            % 读入该表项的坐标值  
            rP=LIS(ls,1);  
            cP=LIS(ls,2);  
            % 根据 Sn 判断该表项‘D’型子孙树是否重要  
            if Sn(1)==1  
                % 每次判断都是读入 Sn 的首位数，判断后立即删除这一位数  
                Sn(1)=\[\];  
                % 生成该表项的孩子树  
                chO=coef\_DOL(rP,cP,'O');  
                % 分别判断每个孩子的重要性  
                for r=1:4  
                    % 读入孩子的坐标值  
                    rO=chO(r,1);  
                    cO=chO(r,2);  
                    % 判断该孩子的重要性  
                    if Sn(1)==1  
                        Sn(1)=\[\];  
                        % 判断该孩子的正负符号  
                        if Sn(1)==1  
                            Sn(1)=\[\];  
                            % 生成该孩子的系数值  
                            DecodeMat(rO,cO)=1.5\*2^N;  
                        else  
                            Sn(1)=\[\];  
                            DecodeMat(rO,cO)=-1.5\*2^N;  
                        end  
                        % 将该孩子添加到重要系数列表 LSP  
                        LSP=\[LSP;chO(r,:)\];  
                    else  
                        % 如果不重要，则这个孩子的系数值为 0  
                        DecodeMat(rO,cO)=0;  
                        Sn(1)=\[\];  
                        % 将该孩子添加到不重要系数列表 LIP  
                        % 本级阈值下不重要的系数在下一级解码中可能是重要的  
                        LIP=\[LIP;chO(r,:)\];  
                    end  
                end  
                % 生成该表项的非直系子孙树  
                chL=coef\_DOL(rP,cP,'L');  
                if ~isempty(chL)  
                    % 如果‘L’型树非空，则将该表项添加到列表 LIS 的尾端等待扫描  
                    LIS=\[LIS;LIS(ls,:)\];  
                    % 表项类型更改为‘L’型  
                    LisFlag=\[LisFlag,'L'\];  
                    % 至此，该表项的‘D’型LIS扫描结束，在LIS中删除该项及其类型符  
                    LIS(ls,:)=\[\];  
                    LisFlag(ls)=\[\];  
                else  
                    % 如果‘L’型树为空集  
                    % 则该表项的‘D’型LIS扫描结束，在LIS中删除该项及其类型符  
                    LIS(ls,:)=\[\];  
                    LisFlag(ls)=\[\];  
                end  
            else  
                % 如果该表项的‘D’型子孙树不重要，则先删除读入的 Sn 值  
                Sn(1)=\[\];  
                % 然后将指针指向下一个 LIS 表项  
                ls=ls+1;  
            end  
            % 更新当前 LIS 的表长，转入下一表项的扫描  
            rlis=size(LIS,1);  
        case 'L'  
            % 对‘L’类表项，不需判断孩子的重要性  
            % 读入该表项的坐标值  
            rP=LIS(ls,1);  
            cP=LIS(ls,2);  
            % 判断该表项的‘L’型子孙树是否重要  
            if Sn(1)==1  
                % 如果该子孙树重要  
                Sn(1)=\[\];  
                % 则生成该表项的孩子树  
                chO=coef\_DOL(rP,cP,'O');  
                % 将该‘L’类表项从 LIS 中删除  
                LIS(ls,:)=\[\];  
                LisFlag(ls)=\[\];  
                % 将表项的四个孩子添加到 LIS 的结尾，标记为‘D’类表项  
                LIS=\[LIS;chO(1:4,:)\];  
                LisFlag(end+1:end+4)='D';  
            else  
                % 如果该表项的‘L’型子孙树不重要，则先删除读入的 Sn 值  
                Sn(1)=\[\];  
                % 然后将指针指向下一个 LIS 表项  
                ls=ls+1;  
            end  
            % 更新当前 LIS 的表长，转入下一表项的扫描  
            rlis=size(LIS,1);  
    end  
end  
% 对 LIS 的扫描结束，将本级阈值的指数减1，准备进入下一级编码  
N=N-1;

4、精细扫描解码程序

function DecodeMat=decRefine(DecodeMat,Rn,N,LSP\_Old)  
% 函数 DECREFINE()为精细解码程序，对上一级解码产生的重要系数列表LSP\_Old，根据输入的  
% 精细位流 Rn 提高重要系数值的重构精度  
% 输入参数：DecodeMat —— 经排序扫描解码后的重构系数矩阵  
%          Rn —— 精细扫描输出位流  
%          N —— 本级解码阈值的指数  
%          LSP\_Old —— 上一级解码产生的重要系数列表  
% 输出参数：DecodeMat —— 提高重要系数精度后的重构矩阵

rlsp=size(LSP\_Old,1);  
if ~isempty(LSP\_Old)  
    for r=1:rlsp  
        dMat=DecodeMat(LSP\_Old(r,1),LSP\_Old(r,2));  
% 首先读取重构矩阵中重要系数的值 dMat  
        rMat=abs(dMat)+(-1)^(1+Rn(1))\*2^(N-1);  
% 对 dMat 的绝对值，如果 Rn = 1，则加上2^(N-1)，否则减去2^(N-1)，结果存入 rMat  
        if dMat<=0  
            rMat=-rMat;  
        end  
% 如果 dMat 为负数，则 rMat 也转为负数  
        Rn(1)=\[\];  
% 消去所读取的Rn信息  
        DecodeMat(LSP\_Old(r,1),LSP\_Old(r,2))=rMat;  
% 将提高了精度的重要系数值返回重构矩阵中  
    end  
end

OK！SPIHT算法的编码、解码程序到这就基本完成啦，下一篇文章我们就用这些程序来做个实例演示。

2、LIP扫描解码程序

function \[DecodeMat,Sn,LSP,LIP\]=lip\_decode(DecodeMat,Sn,N,LSP,LIP)  
% 函数 LIP\_DECODE() 根据排序位流 Sn，更新列表LIP、LSP和重构系数矩阵 DecodeMat  
% 输入参数：DecodeMat —— 上一级解码后生成的重构系数矩阵  
%          Sn —— 本级解码排序位流  
%          N —— 本级解码阈值的指数  
%          LSP —— 上一级解码生成的重要系数列表  
%          LIP —— 上一级解码生成的不重要系数列表  
% 输出参数：DecodeMat —— 本级LIP扫描后更新的重构系数矩阵  
%          Sn —— 对LIP列表扫描后更新的排序位流  
%          LSP —— 对LIP列表扫描后更新的重要系数列表  
%          LIP —— 本级LIP扫描后更新的不重要系数列表

rlip=size(LIP,1);  
% r 是指向 LIP 当前读入表项位置的指针  
r=1;  
% 解码路径与编码路径基本一致，不过解码是根据排序位流 Sn 来判断系数是否重要  
while r<=rlip  
% 读入当前表项的坐标值  
    rN=LIP(r,1);  
    cN=LIP(r,2);  
% 根据 Sn 判断该表项是否重要  
% 根据 Sn 的生成原理，每次判断都是读入 Sn 的首位数，判断后立即删除这一位数  
    if Sn(1)==1  
        % 若Sn(1)=1，则表示当前表项是重要的  
        Sn(1)=\[\];  
        % 读入后即删除该位 Sn 数据，使 Sn(2)变为 Sn(1)，进入下一次判断  
        % 这时的 Sn(1) 是正负符号位  
        if Sn(1)==1  
            % Sn(1)=1，则相应的系数为正数，其值为本级解码阈值的1.5倍  
            DecodeMat(rN,cN)=1.5\*2^N;  
            Sn(1)=\[\];  
        else  
            % Sn(1)=0，则相应的系数为负数  
            DecodeMat(rN,cN)=-1.5\*2^N;  
            Sn(1)=\[\];  
        end  
        % 将该表项添加到重要系数列表 LSP  
        LSP=\[LSP;LIP(r,:)\];  
        % 将该表项从 LIP 中删除  
        LIP(r,:)=\[\];  
    else  
        % 若不重要，则相应的系数值为 0  
        DecodeMat(rN,cN)=0;  
        Sn(1)=\[\];  
        % 将指针指向下一个表项  
        r=r+1;  
    end  
    % 判断当前 LIP 的表长  
    rlip=size(LIP,1);  
end

3、LIS扫描解码程序

function \[LSP,LIP,LIS,LisFlag,DecodeMat,N\]=lis\_decode(DecodeMat,N,Sn,LSP,LIP,LIS,LisFlag)  
% 函数 LIS\_DECODE() 根据排序位流 Sn，更新列表LIP、LSP和重构系数矩阵 DecodeMat  
% 输入参数：DecodeMat —— 本级 LIP 扫描后更新的重构系数矩阵  
%          N —— 本级解码阈值的指数  
%          Sn —— 经本级 LIP 扫描截取后的排序位流  
%          LSP —— 经本级 LIP 扫描后更新的重要系数列表  
%          LIP —— 经本级 LIP 扫描后更新的不重要系数列表  
%          LIS —— 上一级编码生成的不重要子集列表  
%          LisFlag —— 上一级编码生成的不重要子集表项类型列表  
% 输出参数：LSP —— 本级 LIS 扫描后更新的重要系数列表  
%          LIP —— 经本级 LIS 扫描处理后更新的不重要系数列表  
%          LIS —— 本级 LIS 扫描后更新的不重要子集列表  
%          LisFlag —— 本级 LIS 扫描后更新的不重要子集表项类型列表  
%          DecodeMat —— 本级 LIS 扫描后更新的重构系数矩阵  
%          N —— 下一级解码阈值的指数

现在我们分别用一幅16\*16的Lena局部图像和一幅64\*64的Girl局部图像来查看SPIHT算法的编码过程和解码效果。

1、编码过程演示  
    这里我们使用的源图像是一幅16\*16的Lena眼睛局部图像，其源图像和3级小波分解如下图所示：

图1

3级分解的小波系数矩阵为：

利用SPIHT编码算法对图像进行3级编码和3级解码后，输出如下结果：

N =     7            % 初始阈值指数  
T =   128            % 初始阈值

ini\_LSP =     \[\]       % 初始重要系数列表

ini\_LIP =                % 初始不重要系数列表  
     1     1  
     1     2  
     2     1  
     2     2  
     1     3  
     1     4  
     2     3  
     2     4  
     3     1  
     3     2  
     4     1  
     4     2  
     3     3  
     3     4  
     4     3  
     4     4

ini\_LIS =                % 初始不重要子集列表  
     1     3  
     1     4  
     2     3  
     2     4  
     3     1  
     3     2  
     4     1  
     4     2  
     3     3  
     3     4  
     4     3  
     4     4

ini\_LisFlag =    DDDDDDDDDDDD        %初始不重要子集列表中各表项的类型

codeLevel =     1            % 第一级编码  
Threshold =   128           % 第一级编码的阈值  
SnL =                           % 第一级编码输出的排序扫描位流和精细扫描位流  
0  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  7  
RnL =  
7

codeLevel =     2            % 第二级编码  
Threshold =    64           % 第二级编码的阈值  
SnL =                           % 第二级编码输出的排序扫描位流和精细扫描位流  
0  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  7  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  7  
RnL =  
7  0  7

codeLevel =     3            % 第三级编码  
Threshold =    32            % 第三级编码的阈值  
SnL =                            % 第三级编码输出的排序扫描位流和精细扫描位流  
0  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  7  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  7  1  1  0  0  0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0  1  0  0  0  1  0  0  0  0  1  0  0  0  1  0  1  0  0  1  0  0  0  0  0  7  
RnL =  
7  0  7  0  0  1  7

RecImMat =            % 3级解码后重构的小波系数矩阵

可以看到，3级解码后重构的小波系数矩阵大部分区域都是未被扫描赋值的随机数，所以要进行更高级的编码和解码，才能获得高质量的重构图像。  
    对源图像再进行7级编码和7级解码，我们看一下图像重构效果如何：

RecImMat =            % 7级解码后重构的小波系数矩阵

可见，7级解码后系数矩阵大部分元素都被扫描赋值了，只剩下一个2\*2的区域低于解码阈值而未被扫描，但由于初始的重构矩阵是一个元素最大值小于最高级解码阈值的二分之一的随机矩阵，这就保证了重构图像尽可能平滑地恢复，而不会出现区域色块。  
     
2、运行时间  
    下面，我们通过一幅64\*64的Girl图像来看看程序的运行时间和重构效果。  
    图像分解级数：3  
    SPIHT编码级数：7  
    SPIHT解码级数：7

（1）主程序代码如下：

function RecIm=spiht(Im,imDim,codeDim,decodeDim)

global rMat cMat  
strtime=cputime;  
\[rMat,cMat\]=size(Im);  
% ----- Image Wavedec ----- %  
DecIm=mywavedec2(Im,imDim);

% ----- SPIHT Coding ----- %  
\[T,SnList,RnList,ini\_LSP,ini\_LIP,ini\_LIS,ini\_LisFlag\]=spihtcoding(DecIm,imDim,codeDim);

% ----- SPIHT Decoding ----- %  
DecodeMat=spihtdecoding(T,SnList,RnList,ini\_LSP,ini\_LIP,ini\_LIS,ini\_LisFlag,decodeDim);  
% ----- Waverec -----%  
m=2^imDim;  
DecImLL=DecIm(1:rMat/m,1:cMat/m)  
DecodeLL=DecodeMat(1:rMat/m,1:cMat/m)  
% 显示小波分解系数矩阵和重构系数矩阵在LL\_N分解级的数据，比较重构效果  
RecIm=mywaverec2(DecodeMat,imDim,decodeDim);  
Runtime=cputime-strtime

（2）小波分解系数矩阵和重构系数矩阵的LL\_N分解级的数据，比较重构效果

DecImLL =  
   67.0625  106.3750  117.5625  130.3125  110.0000   75.6875   82.6875   65.6250  
   61.5000   50.7500   61.9375  127.0000  104.6875   67.8125   43.2500   36.8750  
   69.3750  117.4375  121.5625  146.7500  123.0625  123.4375   90.9375   92.1250  
   76.9375  153.5000  151.0625  143.3750  111.2500  150.3125  143.3750  124.8125  
   69.8750  135.8750  117.9375  106.4375   85.8750  133.1250  138.6875  116.5625  
   54.2500   94.8750   99.2500  117.9375  118.1250  104.3125  102.0625   88.2500  
   56.2500   75.5625  126.6875  129.4375  119.9375   89.8750  112.0625   83.2500  
   74.5000   45.7500  101.3125  133.0000  124.8125  116.1875  105.0000   77.0000

DecodeLL =  
    67   107   117   131   111    75    83    65  
    61    51    61   127   105    67    43    37  
    69   117   121   147   123   123    91    93  
    77   153   151   143   111   151   143   125  
    69   135   117   107    85   133   139   117  
    55    95    99   117   119   105   103    89  
    57    75   127   129   119    89   113    83  
    75    45   101   133   125   117   105    77

（3）运行时间  
Runtime =    8.0625

最后，我们再看看对普通的256\*256图像，程序运行效果如何。同样是3级小波分解、7级编码、7级解码。

Runtime = 484.8125

转载于:https://www.cnblogs.com/carrey/archive/2009/12/26/1632813.html

[http_blog.csdn.net_chenyusiyuan_archive_2007_12_07_1923334.aspx]: http://blog.csdn.net/chenyusiyuan/archive/2007/12/07/1923334.aspx
[Link 1]: http://blog.csdn.net/chenyusiyuan/archive/2007/12/09/1925291.aspx