本文共 5369 字，大约阅读时间需要 17 分钟。
根据上一篇文章提及的改进算法步骤，编写了如下主要代码，其他的相关程序，例如childMat()、coef_DOL()、SnOut()、fracnum2bin()等函数的代码，请参阅我之前发布的SPIHT算法过程详解与Matlab实现的文章，主要代码如下：
function RecIm=spiht(Im,imDim,codeDim,decodeDim)
global rMat cMat
% ----- 起始时间 ----- %strtime=cputime;
[rMat,cMat]=size(Im);% ----- 小波分解 ----- %DecIm=mywavedec2(Im,imDim);% ----- M_SPIHT Coding ----- %[T,SnList,RnList]=spihtcoding(DecIm,codeDim);% ----- M_SPIHT Decoding ----- %DecodeMat=spihtdecoding(T,SnList,RnList,decodeDim);% ----- 小波重构 -----%RecIm=mywaverec2(DecodeMat,imDim,decodeDim);
% ----- 运行时间 ----- %runtime=cputime-strtime
function [T,SnList,RnList]=spihtcoding(DecIm,codeDim)% 函数 SPIHTCODING() 是SPIHT算法的编码主程序% 输入参数：DecIm ——小波分解系数矩阵；% imDim ——小波分解层数；% codeDim ——编码级数。% 输出参数：T —— 初始阈值，T=2^N，N=floor(log2(max{|c(i,j)|}))，c(i,j)为小波系数矩阵的元素% SnList —— 排序扫描输出位流
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=[];FlagCoef=zeros(rMat,cMat);FlagDch=zeros(rMat/2,cMat/2);FlagLch=zeros(rMat/2,cMat/2);FlagLch(1,1)=1;scanorder=listorder(rMat/2,cMat/2,1,1);%-------------------------%% ----- Coding Loop ----- %%-------------------------%for d=1:codeDim Sn=[]; Rn=[]; % 对系数C(0,0)单独编码 if FlagCoef(1,1)==1 [biLSP,Np]=fracnum2bin(abs(Mat(1,1)),N); Rn=[Rn,biLSP(Np)]; else [isImt,Sign]=SnOut([1 1],N); if isImt Sn=[Sn,1,Sign]; FlagCoef(1,1)=1; else Sn=[Sn,0]; end end % 扫描 FlagDch for i=1:rMat*cMat/4 rD=scanorder(i,1); cD=scanorder(i,2); if FlagDch(rD,cD)==1 chO=coef_DOL(rD,cD,'O'); row=size(chO,1); for j=1:row rO=chO(j,1); cO=chO(j,2); if FlagCoef(rO,cO)==1 [biLSP,Np]=fracnum2bin(abs(Mat(rO,cO)),N); Rn=[Rn,biLSP(Np)]; else [isImt,Sign]=SnOut([rO,cO],N); if isImt Sn=[Sn,1,Sign]; FlagCoef(rO,cO)=1; else Sn=[Sn,0]; end end end end end % 扫描 FlagLch for i=1:rMat*cMat/4 rL=scanorder(i,1); cL=scanorder(i,2); if FlagDch(rL,cL)==0 chD=coef_DOL(rL,cL,'D'); isImt=SnOut(chD,N); if isImt Sn=[Sn,1]; FlagDch(rL,cL)=1; chO=coef_DOL(rL,cL,'O'); row=size(chO,1); for j=1:row rO=chO(j,1); cO=chO(j,2); [isImt,Sign]=SnOut([rO,cO],N); if isImt Sn=[Sn,1,Sign]; FlagCoef(rO,cO)=1; else Sn=[Sn,0]; end end chL=coef_DOL(rL,cL,'L'); if isempty(chL) FlagLch(rL,cL)=0; else isImt=SnOut(chL,N); if isImt Sn=[Sn,1]; FlagLch(rL,cL)=1; chO=coef_DOL(rL,cL,'O'); row=size(chO,1); for j=1:row rO=chO(j,1); cO=chO(j,2); FlagLch(rO,cO)=1; end else Sn=[Sn,0]; end end else Sn=[Sn,0]; end end end N=N-1; SnList=[SnList,Sn,7]; RnList=[RnList,Rn,7]; % 数字‘7’作为区分符,区分不同编码级的Rn、Sn位流end
function DecodeMat=spihtdecoding(T,SnList,RnList,decodeDim)% 函数 SPIHTDECODING() 是SPIHT算法的解码主程序% 输入参数：T —— 初始阈值，T=2^N，N=floor(log2(max{|c(i,j)|}))，c(i,j)为小波系数矩阵的元素% SnList —— 排序扫描输出位流% RnList —— 精细扫描输出位流% decodeDim —— 解码级数% 输出参数：DecodeMat —— 解码后重构的小波系数矩阵
global rMat cMat% rMat、cMat是Mat的行、列数，作为全局变量，在编码、解码的相关程序中使用%-------------------------------------------%% ----- Decoding Input Initialization ----- %%-------------------------------------------%N=log2(T);% 获取初始阈值的指数-NDecodeMat=2^(N-decodeDim)*rand(rMat,cMat);% 初始化重构矩阵为一个随机矩阵，其元素最大值小于最高级解码阈值的二分之一% 这样就可以保证未被扫描赋值的区域有一定的灰度，避免重构图像出现色块
Sn=[];FlagCoef=zeros(rMat,cMat);FlagDch=zeros(rMat/2,cMat/2);FlagLch=zeros(rMat/2,cMat/2);FlagLch(1,1)=1;scanorder=listorder(rMat/2,cMat/2,1,1);%-------------------------%% ----- Decoding Loop ----- %%-------------------------%for d=1:decodeDim [Sn,SnList]=getflow(SnList); [Rn,RnList]=getflow(RnList); % 对系数C(1,1)单独编码 if FlagCoef(1,1)==1 [DecodeMat,Rn]=decRefine(DecodeMat,Rn,N,1,1); else if Sn(1)==1 Sn(1)=[]; if Sn(1)==1 Sn(1)=[]; DecodeMat(1,1)=1.5*2^N; else Sn(1)=[]; DecodeMat(1,1)=-1.5*2^N; end FlagCoef(1,1)=1; else Sn(1)=[]; DecodeMat(1,1)=0; end end % 扫描 FlagDch for i=1:rMat*cMat/4 rD=scanorder(i,1); cD=scanorder(i,2); if FlagDch(rD,cD)==1 chO=coef_DOL(rD,cD,'O'); row=size(chO,1); for j=1:row rO=chO(j,1); cO=chO(j,2); if FlagCoef(rO,cO)==1 [DecodeMat,Rn]=decRefine(DecodeMat,Rn,N,rO,cO); else 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 FlagCoef(rO,cO)=1; else Sn(1)=[]; DecodeMat(rO,cO)=0; end end end end end % 扫描 FlagLch for i=1:rMat*cMat/4 rL=scanorder(i,1); cL=scanorder(i,2); if FlagDch(rL,cL)==0 if Sn(1)==1 Sn(1)=[]; FlagDch(rL,cL)=1; chO=coef_DOL(rL,cL,'O'); row=size(chO,1); for j=1:row rO=chO(j,1); cO=chO(j,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 FlagCoef(rO,cO)=1; else Sn(1)=[]; DecodeMat(rO,cO)=0; end end chL=coef_DOL(rL,cL,'L'); if isempty(chL) FlagLch(rL,cL)=0; else if Sn(1)==1 Sn(1)=[]; FlagLch(rL,cL)=1; chO=coef_DOL(rL,cL,'O'); row=size(chO,1); for j=1:row rO=chO(j,1); cO=chO(j,2); FlagLch(rO,cO)=1; end else Sn(1)=[]; end end else Sn(1)=[]; end end end N=N-1;end
我使用的是Matlab r2007b，小波分解用的是自己编写的Haar二维小波分解程序。对一幅64*64的Girl灰度图像进行SPIHT编码、解码，程序给运行结果：1、按照我前天的文章《多级树集合分裂（SPIHT）算法的过程详解和Matlab实现（2）数学表述》介绍的SPIHT算法步骤编写的程序来进行编、解码，3级小波分解和重构，7级编码，7级解码，程序运行时间是：runtime = 6.8281
2、按照本文所提到的改进算法编写的上述程序来编、解码，同样是3级小波分解和重构，7级编码，7级解码，程序运行时间是：runtime = 9.5313
编解码速度反而下降了，这里的原因是什么呢？欢迎有兴趣的朋友发表评论或者Email交流！