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

上星期我们讨论了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）进行下一级编码扫描。
已经详细介绍了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 —— 下一级编码阈值的指数
global Mat rMat cMat
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用
% 读入当前 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’型子孙列表
global Mat rMat cMat
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用
[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均为空表
global Mat rMat cMat
% Mat是输入的小波分解系数矩阵，作为全局变量，在编码的相关程序中使用
% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用
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编码）。
转自：http://blog.csdn.net/chenyusiyuan/article/details/1925358