2.4 基于小波变换的分形编码算法
2.4.1 小波理论
小波分析是20世纪80年代中期逐渐发展起来的一个数学分支,它与大家熟悉的傅里叶分析有着惊人的相似,其数学思想都来源于经典的调和分析,它是由傅里叶分析逐步发展起来的。
傅里叶变换从提出到至今已经快两个世纪了(1807年由J. Fourier提出)。它是分析信号频率分量的有力工具,但在分析信号的暂态特性时则有点力不从心。1946年,Gabor提出了信号的时频局部化分析方法,即所谓的Gabor变换。后来演变成加窗傅里叶变换(Windows Fourier Transform),也称为短时傅里叶变换。虽然加窗傅里叶变换在一定程度上克服了傅里叶分析的不足,但由著名的Heisenberg测不准原理可知,无论采用什么样的窗函数,该时间函数的时间窗宽度与其傅里叶变换频率窗宽度的乘积不小于ha,b。所以,当选定一窗函数后,使其频宽对应于某一频段时,其时宽就不可能太窄,因而更高频率的信号就不能精确定位。
小波变换的出现使傅里叶分析中遇到的上述问题迎刃而解。小波变换是空间(时间)和频率的局域变换,因而能有效地从信号中提取信息,通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析(Multi-scale Analysis),从而被誉为“数学显微镜”。它不仅继承和发展了窗口傅里叶变换的局部化思想,同时又克服了其不足之处,是一种较理想的对信号进行局部分析的工具,它特别适合于非平稳信号的分析。
小波分析方法的提出可以追溯到1910年Haar提出的小“波”规则正交基及1938年Littlewood-Paley对傅里叶级数建立的L-P理论,按二进制频率成分分组傅里叶变换的相位变化在本质上不影响函数的形状和大小,其后,Caldern于1975年用其早期发现的再生公式给出抛物线空间H′上的原子分解,这个公式后来成了许多函数分解的出发点,它的离散形式已接近小波展开,只是还无法得到组成一正交系的结论。1981年Stromberg对Haar系进行改造,证明了小波函数的存在性。值得注意的是,1984年法国物理学家Morlet在分析地震的局部性质时,发现传统的变换难达到要求,因此,他引入了小波概念于信号分析中,对信号进行分解。随后理论物理学家Grossman尝试对一个确定函数Ψ(x)进行了伸缩、平移操作,有
并对Morlet信号依此方法进行了展开的可行性研究,这无疑为小波分析的形式开创了先河。
真正的小波热开始于1986年,当时Meyer创造性地构造了具有一定衰减性的光滑函数Ψ,其二进制伸缩与平移形式为
是一种构成在L2(R)上的规范正交基。
继Meyer提出小波变换后,Lemarie和Battle又分别独立给出了具有指数衰减的小波函数。1987年,Mallat巧妙地将计算机视觉领域内的多尺度分析的思想引入到小波分析中,小波函数的改造及信号按小波变换的分解与重构,从而成功地统一了前人提出的具体函数的构造,研究了小波变换的离散化情况,提出了著名的Mallat算法,并将其应用于图像的分解与重构。它不仅为纯粹数学的研究提供了有力的工具,而且为图像处理理论的发展树立了里程碑。
小波变换的定义如下。
1.连续小波变换
若基本小波函数为h(x),伸缩和平移因子分别为a和b,则小波是一个满足条件:∫R h(x)dx=0的函数h通过平移和伸缩而产生的一个函数簇ha,b:
通常称为h的基本小波。L2(R)为全实轴上的希尔伯特空间,L2(R)上的可测平方可积函数为f (x),对于任意的f (x)∈L2(R)的连续小波变换定义为
写成内积形式,即有
它对应于f (x)∈L2(R)有函数簇ha,b(x)上的分解,这一分解必须满足下列可容性条件(Admissibility Condition):
其中,H(ω)是h(x)的傅里叶变换,这样的函数h(x)称为解析小波。由式(2.16)可知,函数h(x)可以描述为一带通滤波器的脉冲响应,因此小波变换式(2.14)可以描述为函数f (x)∈L2(R)通过一带通滤波函数的滤波输出,且具有频率放大作用(放大系数为a)。可以证明:如果h(x)满足式(2.16),则由Wa,b(x)可以恢复f (x),即此小波变换的逆变换存在,且有如下定义:
下面给出连续小波变换的两个重要性质。
(1) Parseval公式
设f, g∈L2(R),Wf(a,b)及Wg(a,b)分别表示f、g关于同一基本小波h的小波变换,则有
式中,
(2) 能量公式
设f∈L2(R),Wf(a,b)是f的小波变换,则有
从式(2.17)、式(2.18)和式(2.19)可以看出,小波变换是一种信息保持的可逆变换,原来信号的信息完全保留在小波变换的系数中,变换只是使原信号的能量重新分配。另外,可根据具体问题中f的特征来选取h,使小波变换能更方便地表示f的特征。
下面研究连续小波在相频空间的局部化格式,先假设h是标准的双窗函数,记为
注意到
因此,ha,b的正频窗口中心为(b,)(当a>0时),此外,Δha,b=a · Δh,。
由上可知,尽管由ha,b确定的窗的面积大小固定,但形状各异。对于大的中心频率,窗变窄,即尺寸参数a较小,分辨率在时域或空间域较低,在频域较高;对于小的中心频率,窗变宽,即尺度参数a增大,分辨率在时域或空间域增加,在频域减小。也就是说,在高频端,小波变换在时间上较快,而在低频端,小波变换在频率上较快。这相当于参数a的变化不仅改变连续小波的频谱结构,而且也改变其窗口的大小与形状。
2.离散小波变换
由式(2.15)可通过对其伸缩标度因子a和平移标度因子b的取样而离散化。如果对a和b按如下规律取样:
其中a0>1,b0∈R,m,n∈z2,则由式(2.14)有
这样离散小波变换可定义为
因此,离散小波变换式(2.21)也是一种时频分析,它从集中在某个敬意上的基本函数开始,以规定步长向左或向右移动基本波形,并用标度因子a0来加以扩张或压缩以构造其函数系,一系列小波由此而生,这就是“小波”一词的由来。这里m和n分别称为频率范围指数和时间步长变化指数。由于hm,n正比于,故高频时(对应于小的m值)hm,n高度集中,反之亦然;步长的变化则与n成正比。
为了从离散小波变换式(2.21)重构函数信号f(x),算子
DWm,n:L2(R)→l2(z2)
必须是一有界可逆算子,即对于某个A>0,B<∞,
对一切f (x)∈L2(R)成立。其中 2l(z2)表示二元平方可和序列hm,n(x)向量空间;范数。
3.二进制小波变换
当取a0=2, b0=1及m=j时,将其代入式(2.20),可得到在L2(R)上,该小波正交基为如下形式的函数簇:
Ψj,n(x)={2-j/2Ψ(2-jx-n)(j,n)∈z2}
代入式(2.21),则得到非常重要的离散正交二进制小波变换:
令,函数f (x)可以由它的二进制小波变换重建:
定义2.4 空间L2(R)中的一列闭子空间{Vj}j∈z称为L2(R)的一个多分辨率分析或逼近,如果下面诸条件满足。
① 单调性:对于任意j∈z,Vj⊂Vj-1。
② 逼近性:,。
③ 伸缩性:u(x)∈Vj⇔u(2x)∈Vj1-。
④ Riesz基:存在φ∈V0使{φ(x-k)k∈z}构成V0的标准正交基。
⑤Vj=wj+1⊕wj+2⊕wj+3…。
条件③表明,空间列{Vj}中任一空间Vi的基经过简单的伸缩变换而得到。因此,若{φ(x-k)ㄧk∈z}是V0的标准正交基,则对于任意j∈z,函数基{φj,n(x)=2-j/2φ(2-jx-n)ㄧn∈z}构成Vj的标准正交基。函数φ(x)称为尺度函数,具有低通滤波的作用,且满足:
其中,
式(2.22)称为双尺度方程。
在Vi的正交补空间wj中,函数基{Ψj,n(x)=2-j/2(Ψ2-jx-n)n∈z}构成wj的标准正交基。函数 (x)Ψ 称为小波函数,具有高通滤波的作用。其中g(n)=(-1)nh(1-n)。在信号处理中,g(n)和h(n)称为正交镜像滤波器。
2.4.2 基于小波变换的图像分形编码算法
采用小波变换出于两个方面的考虑,其一,图像的细节变化和少许遮掩只影响图像中高频部分的变化,这样图像的低频部分就在有细节变化的情况下仍然比较稳定,因此可以考虑仅对图像的低频部分进行编码,而利用小波变换可以达到这个目的;其二,对整幅图像进行编码计算量较大,而如果只利用小波变换的低频图像的话,数据量就减少到了原来的四分之一(仅做一次变换),可以较大的提高运算速度。具体实现时是对原人脸图像做一次小波变换,对得到的低频图像做进一步处理。
基于小波多分辨分析的图像分解算法的基本思想是,在选定合适的正交小波基的基础上,对图像进行一阶小波分解,得到对应于图像的一个低频子图和三个高频子图;图像信号的主要能量集中在低频区域,它反映图像的平均亮度,所以只须对这一部分进行四叉树分形编码。实验结果表明,该方法能够保证较高的信噪比,提高了人脸图像分形编码质量。图2.8给出了对原图像进行小波分解的实例。
图2.8 对原图像进行小波分解
基本自动分形图像编码(Basic Automatic Fractal Image Coding, BAFIC)方法中,图像分块大小固定,不能很好反映分形特征。若用较大固定块(如8×8个像素),则图像中复杂度较高区域(如眼睛等)得不到很好的匹配,从而解码图像质量将受到影响。若将分块尺寸减小(如4×4),则能得到较好的匹配。但图像中很多较为平滑区域(如背景部分),用较大分块就能得到很好匹配,再把区域分块,势必增加编码所需映射数目,从而降低压缩率并增加计算量。为解决这一矛盾,就必须寻找更合理的图像分割方法。采用四叉树方法对一阶小波分解的低频子图进行分割,从而较好地解决了上述矛盾。二值图像四叉树表示,是基于数据区域一致性判别准则和空间递归分解原理建立起来的树形数据结构。分解过程可用四叉树表示,根结点表示整幅图像,叶结点表示黑或白四分区,非叶结点表示灰四分区。图2.9是图像分解及其四叉树表示。
图2.9 图像分解及其四叉树表示
在BAFIC方法中,Range块在寻找匹配Domain块过程中计算量非常大,大大增加了编码时间,所以有必要对Range块、Domain块进行分类,以便减少与一个Range块相比较的Domain块的数目。在编码之前Domain集中所有块都被分类。编码过程中,再对正待编码的Range块分类。只有和该Range块在同一类(或相近类) Domain块才和Range块进行匹配比较。
经过上述分析,给出新的分类方法[2.22]如下:先把一个图像块分成NW、NE、SW、SE四个等大小的四分块,顺序标号为1、2、3、4。对每个四分块计算其像素的平均值和方差,即若第i(i=1, 2, 3, 4)个四分块的像素值为R1′,R2′,…,Rn′,则计算
然后对图像块进行旋转和翻转,以使Ai的排列顺序为下面三种之一。
主类1:A1≥A2≥A3≥A4。
主类2:A1≥A2≥A4≥A3。
主类3:A1≥A4≥A2≥A3。
它们对应的亮度如图2.10所示。一旦固定了图像子块的方向,根据Vi的大小顺序,每个主类中还包含24个子类。这样一共有72类图像块。
图2.10 主类亮度顺序图示
编码步骤如下:
① 原始标准人脸图像A,分割成不相交叠的个子块Ri作为初始的分割,Ri的左上角位于(K,L)的位置上,用取块符A表示。
② 对92×112像素的ORL中的人脸图像,设定最小深度、最大深度分别是3和5。
③ 对正待编码的Range块分类,只有和该Range块在同一类(或相近类)Domain块才和Range块进行匹配比较。
④ 搜索其相应的父块,位置在π(K,L)上,采用“反射-旋转”操作,和平均抽样操作搜索的相应父块为LP(K,L)AvA。设子块与变换后的父块的灰度分别为1,aa2,…,am b1,b2,…,bm。子块与变换后的父块之间的误差为,
令∂dp ∂CK,L=0, ∂dp ∂hK ,L=0,求出, ,按下式计算出。
⑤ 如果,则记下(K,L)、π(K,L)、P(K,L)、、,就完成了对Ri的编码,否则,转到④。
⑥ 如果,则把子块Ri再划分为四块,父块也相应划分为四块,此时取t=t +1,子块仍用tA表示。对父块的变换相应为LP(K,L)AvA,继续采用①,此时,
这里的bi′、ai′都是在四叉树划分时块中的灰度值,并按式(2.25)求出。
⑦ 重复③和④,直到所有的大小子块与相适应的父块之间自相似变换误差都小于ε为止。
编码文件是由参数和一些有效数据构成。参数包括:rsm误差门限、四叉树最小和最大深度级、Domain集类型及步长参数、允许最大灰度比例系数、待编码图像尺寸、与Range块相比较的Domain类的数目。有效数据如下:图像的最终四叉树结构、每个Range块的灰度比例系数和灰度偏移系数及它的最佳匹配Domain块、映射Domain块到Range块的方向信息。
为了提高压缩率,使用存储分配方案为:
① 每个不同树深度的图像块分配1 bit来表明是否还要继续细分。
② 灰度比例系数5 bit,灰度偏移系数7 bit。
③ 每个Domain块,要记录其位置。但是,当灰度比例系数为0时,就和Domain块无关,因此也就无须存储Domain及方向信息。
④ Domain-Range映射可能的方向有8种,要用3 bit来存储这种信息。
解码过程是在任意一幅初始图像上进行迭代,迭代后的解码图像由若干迭代变换所构成的分形图像(对应于编码时的Range块,这些Range块由四叉树分割产生)组成。分形有一个重要特性就是包含任意比例的细节,所以解码图像与分辨率无关。