多级树集合分裂（SPIHT）算法的过程详解和Matlab实现（3）编码过程——主程序

上一篇文章已经详细介绍了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,:);