多级树集合分裂SPIHT算法的过程详解与Matlab实现.docx
- 文档编号:5936443
- 上传时间:2023-01-02
- 格式:DOCX
- 页数:22
- 大小:501.79KB
多级树集合分裂SPIHT算法的过程详解与Matlab实现.docx
《多级树集合分裂SPIHT算法的过程详解与Matlab实现.docx》由会员分享,可在线阅读,更多相关《多级树集合分裂SPIHT算法的过程详解与Matlab实现.docx(22页珍藏版)》请在冰豆网上搜索。
多级树集合分裂SPIHT算法的过程详解与Matlab实现
多级树集合分裂(SPIHT)算法的过程详解与Matlab实现
上星期我们讨论了EZW算法,很高兴收到了一些朋友的email,对算法进行探讨、交流。
这也是我开这个博客的源动力之一,学习就应该开诚布公、交流互助,在探讨中加深对所学知识的理解和掌握。
在弄懂了EZW算法原理并用Matlab实现后,我继续学习EZW的改进算法。
至今有一周的时间没更新博客、写新文章了,其实就是把时间用在EZW的一个改进算法——多级树集合分裂(SetPartitioninginHierarchicalTrees,SPIHT)算法的研究和Matlab实现。
由于EZW是SPIHT的基础,所以在EZW算法的Matlab代码的基础上,我很快就完成了SPIHT的代码编写,但最痛苦的是一开始没吃透算法原理,程序在初始分集上出了错,调试了两天找不出根本问题,昨天从头再看一次算法原理,才发现问题所在……呵呵,小小的粗心就耽搁了我两三天的时间和精力!
问题解决后,就编写程序注释了,上次EZW算法的代码都没写注释,让大家看着辛苦,不好意思哦!
好,接下来就开始讨论SPIHT算法的原理,然后给出具体的Matlab代码。
一、SPIHT算法与EZW算法
EZW算法是一种基于零树的嵌入式图象编码算法,虽然在小波变换系数中,零树是一个比较有效的表示不重要系数的数据结构,但是,在小波系数中还存在这样的树结构,它的树根是重要的,除树根以外的其它结点是不重要的。
对这样的系数结构,零树就不是一种很有效的表示方法。
A.Said和W.A.Pearlman根据Shapiro零树编码算法(EZW)的基本思想,提出了一种新的且性能更优的实现方法,即基于多级树集合分裂排序(SetPartitioninginHierarchicalTrees,SPIHT)的编码算法。
它采用了空间方向树(SOT:
spatialorientationtree)、全体子孙集合D(i,j)和非直系子孙集合L(i,j)的概念以更有效地表示上述特征的系数结构,从而提高了编码效率。
SPIHT算法能够生成一个嵌入位流(embeddedbitstream),使接收的位流在任意点中断时,都可解压和重构图像,具有良好的渐进传输特性;算法的初始化过程、细化过程类似于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型两种
globalMatrMatcMat
%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)为小波系数矩阵的元素
%--------------------------------------%
%-----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=[];
%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型。
%------------------------------------------------%
%-----CodingInputInitialization------%
%------------------------------------------------%
LSP=ini_LSP;
LIP=ini_LIP;
LIS=ini_LIS;
LisFlag=ini_LisFlag;
%将待输出的各项列表存入相应的编码工作列表
%--------------------------------%
%-----CodingLoop------%
%--------------------------------%
ford=1:
codeDim
%------------------------------------------%
%-----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
%各项输出参数均作为下一编码级的输入
%-------------------------------------%
%-----RefinementPass-----%
%-------------------------------------%
Rn=refinement(N+1,LSP_Old);
%精细扫描是在当前阈值T=2^N下,扫描上一编码级产生的LSP,故输入为(N+1,LSP_Old),
%输出为精细位流Rn
%-----------------------------------%
%-----OutputDataflow-----%
%-----------------------------------%
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的程序
functionlp=coef_H(imDim)
%函数COEF_H()根据矩阵的行列数rMat、cMat和小波分解层数imDim来计算树根坐标集H
%输入参数:
imDim——小波分解层数,也可记作N
%输出参数:
lp——rMat*cMat矩阵经N层分解后,LL_N,LH_N,HL_N,HH_N所有系数的坐标集合
globalrMatcMat
%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列表,以下是该函数的程序代码:
functionlsorder=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编码的精细扫描程序,其中包括一个能够将带小数的十进制数转换为二进制表示的函数,这个转换函数可以实现任意精度的二进制转换,特别是将小数部分转换为二进制表示。
希望对有需要的朋友有所帮助。
下一篇文章将给出SPIHT的解码程序。
请关注后续文章,欢迎Email联系交流。
4、精细扫描程序
functionRn=refinement(N,LSP_Old)
%函数REFINEMENT()为精细编码程序,对上一级编码产生的重要系数列表LSP_Old,读取每个
%表项相应小波系数绝对值的二进制表示,输出其中第N个重要的位,即相应于2^N处的码数
%输入参数:
N——本级编码阈值的指数
% LSP_Old——上一级编码产生的重要系数列表
%输出参数:
Rn——精细扫描输出位流
globalMat
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 多级 集合 分裂 SPIHT 算法 过程 详解 Matlab 实现