第14章 小波分析用于图像压缩

图像压缩的基本理论起源于20世纪40年代末香农(Shannon)的信息理论。香农定理表明，在不产生任何失真的前提下，通过合理的编码，对于每一个信源符号分配不等长的码字可以任意接近于信源的熵。在这个前提下出现了几种不同的无失真信源编码方法，如huffman编码、算术编码、词典编码等，这些方法可以应用于一幅数字图像，能获得一定的码率压缩。但是无失真压缩是很有限的，对较复杂的自然图像，压缩率一般不超过两倍。
本章从实践角度出发，首先介绍图像压缩的基本方法以及在工程实践中的压缩标准；其次，介绍有关MATLAB小波工具箱中的压缩函数，并且将不惜篇幅地列举压缩的算例供读者参考；最后，介绍利用小波分析进行图像压缩的综合实例。
学习目标：
(1)了解小波压缩的原理和方法
(2)熟练掌握MATLAB中小波压缩函数
(3)熟练掌握二维小波工具箱对图像压缩
14.1 图像压缩介绍
与图像去噪相似，压缩领域中由于小波的特殊优点，应用其进行压缩也受到了许多学者的重视，并获得了非常好的效果。
14.1.1 数据冗余
图像压缩所要解决的问题是尽可能减少表示数字图像需要的数据量。减少数据量的基本原理是除去其中的多余数据。以数学的观点来看，这一过程实际上就是将二维像素阵列变换为一个在统计上没有关联的数据集合。
这种变换在图像存储或者传输之前进行，在以后的某个时候，再对压缩图像进行解压来重构原图像或原图像的近似图像。
数据是用来表示信息的，如果不同的方法为表示给定量的信息使用了不同的数据量，那么使用较多数据量的方法中，有些数据必然是代表了无用的信息，或者是重复地表示了其他数据已经表示了的信息，这就是数据冗余的概念，它是数据压缩中的关键。
在数字图像压缩中，有3种基本的数据冗余：
(1)编码冗余；
(2)像素间冗余；
(3)心理视觉冗余。
如果能够减少或者消除其中的一种或者多种冗余，就能取得数据压缩的效果。
1.编码冗余
对图像编码需要建立码本以表达图像数据。这里的码本是指用来表达一定量的信息或者一组事件所需的一系列的符号(如字母、数字等)。其中对每个信息或时间所赋的码符号序列称为码字，而每个码字的符号个数称为码字长度。
我们以一幅256×256的8位灰度图为例：其每个像素的灰度值在0～255之间，设k表示0～255之间的灰度值，Nk表示灰度值为k的像素总个数，N为图像的总像素数，P(k)表示像素灰度值为k的像素出现的概率：P (k)=N (k)/N，k∈[0，255]。
设用来表示灰度值k的每个数值的比特数是L(k)，那么为表示每个像素所需的平均比特数是：
最简单的二元码本成为自然码。对每个信息或事件所赋的码是从2m
 个m bit的二元码中选出来的一个。如果用自然码表示一幅图像的灰度值，则由上式得出平均码长为m。
根据上式，如果用较小的比特数表示出现概率较大的灰度级，而用较多的比特数表示出现概率较小的灰度级，就能达到数据压缩的效果。这种压缩常称为变长码。如果编码所用的码本不能使式达到最小，则说明存在编码冗余。
一般来说，如果编码时没有充分利用编码对象的概率特性就会产生编码冗余。
2.像素间冗余
在图像中，像素间的冗余通常称为空间冗余或者几何冗余。各像素之间的值可以比较方便地由其临近的像素表示出来，每个独立的像素所携带的信息相对较少。换句话说，单个像素对图像的视觉贡献有很多是冗余的，因为常能用基于其临近像素的值来推断。
为了减少图像中的像素冗余，需要将常用的2D像素矩阵表达式转换为某种更有效(但可能不直观)的表达形式。这种减少像素间冗余的转换常称为映射。如果原始图像元素能从转换后的数据集合重建出来，则这种映射称为可反转的，否则就称为不可反转的。
3.心理视觉冗余
眼睛所感受到的图像区域亮度不仅仅与区域的反射光有关，例如马赫带效应，在灰度值为常数的区域也能感觉到灰度值的变化。这种现象的产生是由于眼睛并不是对所有视觉信息有相同的灵敏度。有些信息在通常的视觉过程中与另外一些信息相比来说并不是那么重要，这些信息可以认为是心理视觉冗余的，去除这些信息并不会明显地降低所感受到的图像质量。许多称为“第二代编码技术”的方法就是基于这个原理的。
心理视觉冗余的存在是与人观察图像的方式有关的。人在观察图像时主要是寻找某些比较明显的特征，而不是定量的分析图像中每个像素的亮度，或至少不是对每个像素等同的分析。人通过在脑子里分析这些特征并与先验知识结合以完成对图像的解释过程。由于每个人所具有的先验知识不同，对同一幅图的心理视觉冗余也是因人而异。
心理视觉冗余从本质上说与前两种冗余不同，它是与实在的视觉信息联系着的。只有在这些信息对正常的视觉来说并不是必不可少的时候才可能被去除。因为去除心理的视觉冗余数据能导致定量信息有损失，所以量化是不可逆转操作，它用于数据压缩会导致有损压缩。根据心理视觉冗余的特点，可以采取一些有效的措施来压缩数据量，电视中的隔行扫描就是一个例子。
14.1.2 变换编码
由上一小节所描述的图像数据存在很多冗余，但是如果直接在空域直接处理比较困难，通常通过变换将图像数据从空域变换到变换域，对变换域进行处理。
对变换的要求通常是要将强相关的空间域的数据映射到完全不相关的、能量分布集中的变换系数矩阵，占少数的大的变换系数代表了图像的最主要的能量成分，占多数的小的变换系数表示了一些不重要的细节分量。同时变换要能够容易的与人们关于人类视觉系统的知识相结合，以有效地取出视觉冗余，尽可能地保留重要视觉信息。
具备最理想的去相关和能量集中紧致性最强的变换是KL变换，但是KL变换的复杂性使得其应用不现实。由于DCT变换能够很好地逼近KL变换，同时运算比较简单，因此DCT广泛地应用于图像压缩，一直成为变换编码的主要工具，例如JPEG标准就是采用的DCT变换。
20世纪80年代后期，小波变换的发展提供了一种新的变换工具，小波变换能够适应现代图像压缩所需要的如多分辨率、多层质量控制等要求，在较大压缩比下，小波图像压缩质量明显好于DCT变换，因此在新的JPEG 2000标准中也采用了小波变换作为核心算法。
14.1.3 图像压缩模型
现在我们来介绍一个通用的图像压缩编码模型，如图14-1所示。
图14-1 通用压缩编码模型
图像压缩目的是减少或消除图像中的编码冗余、像素间冗余及心理视觉冗余。在图像压缩系统中，映射器将输入数据变换以减少像素间冗余。这个操作一般是可以逆转的，它可以直接减少也可以不直接减少图像数据，这与具体编码技术有关。变换编码中的变换算法就是一种映射器，通过某种数学变换将图像数据从一个域转换至另一个域，以使图像从复杂的表示变成一种更简单的或者说是具有更好统计特性的表示，减少数据的相关性。
量化器根据给定的保真度准则减少映射器输出的精确度。这个操作可以减少心理视觉冗余，但是不可反转，所以不可用在无损压缩中。符号编码器产生表达量化器输出的码本，并根据码本映射输出。一般情况下采用变长码来表达映射和量化后的数据。它通过将最短的码赋给最频繁出现的输出值以减少编码冗余。这个操作是可以反转的。
需要指出的是，不是所有的图像压缩编码器都一定包含以上 3个子模块。例如无失真信源编码器就不能有量化器，另外有些压缩技术常把上述物理上可分离的子模块结合起来。
14.1.4 图像压缩技术
对声音、图像、视频等多媒体信息的压缩有两条思路，要么采用成熟的通用数据压缩技术进行压缩，要么根据媒体信息的特性设计新的压缩方法。
GIF是使用通用压缩技术压缩图像信息的最成功的例子之一，它使用LZW压缩算法，可以把原始图形文件以非常小的数据量存储，可以在同一个文件中存储多幅图像从而实现动画效果。
GIF文件中除了经过LZW压缩的像素信息以外，还保存有图像的各种属性信息以及图像所使用的调色板信息等。GIF 精确地保留了原始图像的每一个像素信息，是无损图像压缩的代表。因为GIF的成功，直到现在这种格式仍然广为应用。
根据媒体特性量身定制的压缩方法中，行程编码(Run-Length Encoding，RLE)是最为简单、最容易被想到的一种。大多数计算机中产生的图像(和现实世界的图像例如照片不同)都具有着大面积重复的颜色块，可以用一个颜色值加一个重复次数来表示这一块图像，冗余度由此减小了，这就是RLE方法的基本思路。但是，它不适于用来压缩照片、声音等连续重复信息很少的数据。RLE方法最有代表性的实现有PCX和TIFF图形格式，其最为成功的应用领域是传真机中传真图像的压缩。
20世纪70年代末80年代初，人们逐渐意识到，对大多数灰度或是彩色图像乃至声音文件，没有必要忠实地保留其所有信息，在允许一定的精度损失的情况下，可以实现更为有效的压缩方法。到20世纪80年代末，许多人已经在这一领域取得了不小的成果，设计出了一批压缩效果较好的声音和图像压缩算法。
在此基础上，国际标准化组织(ISO)和CCITT联合组成了两个委员会：静态图像联合专家小组(JPEG)和动态图像联合专家小组(MPEG)。JPEG的压缩目标是静止图像(灰度的和彩色的)，MPEG的目标则是声音和视频。但他们的基本思路是完全一样的，即保留媒体信息中最有规律、最能体现信息主要特征的数据，而略去其他不重要的数据。他们都取得了令人赞叹的成就。
JPEG(Joint Photographic Experts Group)是在国际标准化组织(ISO)领导之下制定静态图像压缩标准的委员会。第一个国际的静态图像压缩标准——ISO 109110-1(JPEG)就是该委员会制定的。
由于 JPEG 所具有的优点，它在短短几年内获得了极大的成功，被广泛应用于互联网和数码相机领域，网站上80%的影像都采用了JPEG压缩标准。
JPEG标准主要涉及连续色调(灰度和彩色)静止图像的压缩编码，共提供了4种工作模式，如下。
•顺序型模式：即按照通常的从左到右、从上到下的顺序对图像进行编码，称为基本系统。
•分层模式：按照分辨率由低到高进行编码。
•渐进模式。
•无失真模式。
其中，JPEG基本系统是该标准的核心，提供了一种简单高效的图像编码方案。JPEG 2000采用改进的压缩技术来提供更高的解像度，其伸缩能力可以为一个文件提供从无损到有损的多种画质和分辨率选择。它放弃了 JPEG 所采用的以离散余弦变换(DCT)算法为主的区块编码方式，而改用以离散小波变换算法为主的多解析编码方式。
在压缩率相同的情况下，JPEG 2000的信噪比将比JPEG平均提高2dB左右。JPEG 2000拥有3种层次的编码形式：彩色静态画面采用的JPEG编码、2值图像采用的JBIG、低压缩率图像采用的JPEGLS等，成为应用于各种图像的通用编码方式。
在编码算法上，JPEG 2000采用离散小波变换(DWT)和Bit Plain算术编码(MQ coder)。JPEG中使用的离散余弦变换算法(DCT)是经典谱分析常采用的工具，它考察整个时域过程的频域特征或整个频域过程的时域特征。因为在处理时图像数据首先被分割为8×8的小块进行处理，所以恢复出来的图像有分块效应(即马赛克效应)，压缩比越大，分块效应越明显。
JPEG 2000中使用的离散小波变换算法(DWT)是现代谱分析工具，对高频成分采用由细到粗渐进的时空域上的间隔取样，所以能够像变焦摄影一样看清远近不同的景物，并放大任意细节，是构造图像多重分辨率的有力工具。
与传统JPEG压缩方式对比，在较高压缩比率下，JPEG 2000的表现更为出色。用JPEG压缩的图像，在压缩比较高的情况下，有明显的马赛克现象(比如出现较大色斑或颜色信息丢失)，但是用JPEG 2000 压缩的图像效果就能得到保证。即使在很高的压缩比下，图像的内容也能够辨别。
同时，JPEG 2000 纠错能力很强，并且能够在非标准的图像尺寸下保证再现较高图像质量。
14.1.5 JPEG 2000压缩算法
目前的 JPEG 静止图像压缩标准，在中高比特率上压缩效果较好，然而在低比特速率的情况下，重构图像存在严重的方块效应，不能很好地适应网络传描图像的需要。尽管目前的JPEG标准有44种操作模式，但是其中的大部分模式是针对不同的应用提出的，不具有通用性，这给交换、传输的压缩图像带来很大的麻烦。
为了弥补目前标准的不足，适应21世纪图像压缩的需要，早在1997年ISO/ITU-T组织下的IECJTCl/SC29/WGl 小组便开始着手制定新的静止图像压缩标准——JPEG 2000。与JPEG不同，JPEG 2000基于小波变换，采用当前最新的嵌入式编码技术，在获得优于目前JPEG标准压缩效果的同时，生成的码流有较强的功能，可应用于多个领域。
1.小波变换
不同于传统的DCT变换，小波变换具有对图像进行多分辨分析和反映图像局部特征的特点。通过对图像进行离散小波变换，得到小波系数图像，而分解的级数视具体情况而定。
小波系数图像由几种子带系数图像组成。这些子带系数图像描述的是图像水平和垂直方向的空间频率特性。不同子带的小波系数反映图像不同空间分辨率的特性。通过多级小波分解，小波系数既能表示图像中局部区域的高频信息(如图像边缘)，也能表示图像中的低频信息(如图像背景)。这样，即使在低比特率的情况下，我们也能保持较多的图像细节(如边缘)。另外，下一级分解得到的系数所表示图像在水平和垂直方向的分辨率只有上一级小波系数所表示的图像的一半。所以，通过对系数图像的不同级数进行解码，就可以得到具有不同空间分辨率(或清晰，或模糊)的图像。
小波变换因其具有的这种优点被JPEG 2000 标准所采用。在编码系统中，对每个图像进行Mallat塔式小波分解。经过大量的测试，JPEG 2000 选用两种小波滤波器：LeGall 5/3滤波器和Daubechies 9/7 滤波器。前者可用于有损或无损图像压缩，后者只能用于有损压缩。
在JPEG 2000 标准中，小波滤波器可以有2 种实现模式：基于卷积的和基于提升机制的。而具体实现时，对图像边缘都要进行周期对称延伸，这样可以防止滤波器对图像边缘操作时产生失真。另外，为了减小变换时所需空间的开销，标准中还应用了基于行的小波变换技术。
2.量化
由于人类视觉系统对图像的分辨率要求有一定的局限，通过适当的量化减小变换系数的精度，可在不影响图像主观质量的前提下达到图像压缩的目的。量化的关键是根据变换后图像的特征、重构图像质量要求等因素设计合理的量化步长。
量化操作是有损的，会产生量化误差。不过一种情况除外，那就是量化步长是1，并且小波系数都是整数，利用可恢复整数 5/3 拍小波滤波器进行小波变换得到的结果就符合这种情况。
在JPEG 2000 标准中，对每一个子带可以有不同的量化步长，但是在一个子带中只有一个量化步长。量化以后，每一个小波系数有2部分来表示：符号和幅值。对量化后的小波系数进行编码。对于无损压缩，量化步长必须是1。
3.熵编码
图像经过变换、量化后，在一定程度上减少了空域和频域上的冗余度，但是这些数据在统计意义上还存在一定的相关性，为此采用熵编码来消除数据间的统计相关。将量化后的子带系数划分成小的矩形单元——码块(Code Block)。
4.位流组织
为了适合图像交换，更好地应用JPEG 2000 压缩码流的功能，JPEG 2000 标准规定了存放压缩位流和解码所需参数的格式，把压缩码流以包为单元进行组织，形成最终的码流。
JPEG 2000 通过采用速率控制方法来计算码流的理想截断点，从而获得给定压缩码率下的最佳重构图像质量。
速率控制使用了PCRD率失真优化算法。率失真优化，即给定整个压缩码流的最大编码速率，找出每个码块压缩位流的适当截断点，在满足的条件下，使重构图像失真最小。
从而使得嵌入式码块编码具有如下特点：生成的压缩位流可根据需要被截断成不同长度的位流子集，将所有码块的截断位流组织起来可重构出一定质量的图像。
14.1.6 JPEG与JPEG 2000的区别
JPEG 2000与传统 JPEG 最大的不同在于它放弃了 JPEG 所采用的以离散余弦变换(DCT)为主的区块编码方式，转而采用以小波变换(DWT)为主的多解析编码方式。
余弦变换是经典的谱分析工具，它考察的是整个时域过程的频域特征或整个频域过程的时域特征，因此对于平稳过程，它有很好的效果，但对于非平稳过程，它却有诸多不足。在JPEG中，离散余弦变换将图像压缩为8×8的小块，然后依次放入文件中，这种算法靠丢弃频率信息实现压缩，因而图像的压缩率越高，频率信息被丢弃的越多。在极端情况下，JPEG图像只保留了反映图像外貌的基本信息，精细的图像细节都损失了。
小波变换是现代谱分析工具，它既能考察局部时域过程的频域特征，又能考察局部频域过程的时域特征，因此即使对于非平稳过程，处理起来也得心应手。它能将图像变换为一系列小波系数，这些系数可以被高效压缩和存储，此外，小波的粗略边缘可以更好地表现图像，因为它消除了DCT压缩普遍具有的方块效应。
现在网络上的JPEG图像下载时是按“块”传输的，因此只能逐行显示，而采用JPEG 2000格式的图像支持渐进传输。JPEG 2000中的渐进传输有2 种：按照分辨率的渐进传输和按照质量的渐进传输。
按照质量的渐进传输就是先传输图像轮廓数据，然后再逐步传输细节数据来不断提高图像质量，而按照分辨率的渐进传输则先传输分辨率较低的图像，后一幅图像在前一幅图像的基础上提高其分辨率。图像的渐进传输使得用户不需要像以前那样等图像全部下载后才决定是否需要，有助于快速地浏览和选择大量图片，从而有效地解决了网络传输的不通顺问题。
JPEG 2000一个极其重要的优点就是感兴趣区域(Region of Interest，ROI)。用户可以任意指定图片上感兴趣的区域，然后在压缩时对这些区域指定压缩质量，或在恢复时指定某些区域的解压缩要求。这是因为子波在空间和频率域上具有局域性，要完全恢复图像中的某个局部，并不需要所有编码都被精确保留，只要对应它的一部分编码没有误差就可以了。
在实际应用中，我们就可以对一幅图像中感兴趣的部分采用低压缩比以获取较好的图像效果，而对其他部分采用高压缩比以节省存储空间。这样就能在保证不丢失重要信息的同时又有效地压缩了数据量，实现了真正的“交互式”压缩。
在相同的压缩质量参数下，JPEG 2000的峰值信噪比要优于JPEG；且在低位率时JPEG压缩会出现明显的方块效应，JPEG 2000 则会出现Gibbs 效应，但JPEG 2000的保真度要优于JPEG。
JPEG 2000 作为JPEG的升级版，其压缩率比JPEG高约30%，同时支持有损和无损压缩。JPEG 2000格式有一个极其重要的特征在于它能实现渐进传输，即先传输图像的轮廓，然后逐步传输数据，不断提高图像质量，让图像由朦胧到清晰显示。
JPEG 2000 旨在创建一个新的图像编码系统，该压缩编码系统的率失真和主观图像质量优于现有的 JPEG 标准，能够提供对图像的低码率的压缩，并且对压缩码流可进行灵活的处理，如随机获取部分压缩码流、图像的渐进传输、感兴趣区域的实现以及压缩码流较强的容错性能等，该标准将与现行JPEG标准兼容。
JPEG 2000和JPEG相比优势明显，且向下兼容，因此可取代传统的JPEG格式。JPEG 2000既可应用于传统的JPEG市场，如扫描仪、数码相机等，又可应用于新兴领域，如网路传输、无线通信等。
JPEG不能在单一码流中提供有损和无损压缩，并且不能支持大于64×64KB的图像压缩。同时，尽管当前的 JPEG 标准具有重新启动间隔的规定，但当碰到比特差错时图像质量将受到严重的损坏。
虽然JPEG 2000 比JPEG 有很多优点，但是由于它的技术还不成熟，还存在稳定性、复杂性和准备不充分等问题，所以目前基于DCT的JPEG图像压缩编码仍然是世界上图像压缩技术标准应用最广泛的标准。
14.1.7 基于DCT的JPEG图像压缩编码
JPEG 文件使用的颜色空间为1982年推荐的电视图像数字化标准 CCIR 601(现为ITU-RB T.601)。在这个色彩空间中，每个分量、每个像素的电平规定为255 级，用8 位代码表示。JPEG只支持YCbCr颜色模式。其中Y 代表亮度，CbCr代表色度。全彩色图像RGB模式转换到YCbCr模式，用下组公式：
其逆变换为：
JPEG是以8×8的块为单位来进行处理的，由于人眼对亮度Y的敏感度比色度CbCr的敏感度大得多，所以采用缩减取样的方式，通常采用YUV422取样，如图14-2所示。即对于16×16的块，Y取4个8×8的块，Cb、Cr各取2个8×8的块。也有YUV411方式，Y取4个8×8的块，Cb、Cr各取1个8×8的块。YUV422取样方式，数据减少1/3。YUV411取样方式，数据减少1/2。缩减取样一般采用图14-3所示方法。
图14-2 YUV422取样示意图
图14-3 压缩取样示意图
在傅立叶级数展开式中，如果被展开的函数是实偶函数，那么，其傅立叶级数中只包含余弦项，在将其离散化由此可导出余弦变换，或称之为离散余弦变换(Discrete Cosine Transform，DCT)。
二维离散余弦正变换公式为：
式中：
二维离散余弦逆变换公式为：
式中：
JPEG采用的是8×8大小的子块的二维离散余弦变换。在编码器的输入端，把原始图像顺序地分割成一系列8×8的子块，子块的数值在−128～127之间。采用余弦变换获得64个变换系数。
变换公式，如式(14.5)所示：
式中：x，y，u，v=0，1，"，7。
1.DCT系数的量化
量化是对经过DCT变换后的频率系数进行量化，其目的是减小非“0”系数的幅度以及增加“0”值系数的数目，它是图像质量下降的最主要原因。
对于基于DCT的JPEG图像压缩编码算法使用如图14-4所示的均匀量化器进行量化，量化步距是按照系数所在的位置和每种颜色分量的色调值来确定的。因为人眼对亮度图像比对色差图像更敏感，因此使用了表14-1所示的量化表。
图14-4 均匀量化器
表14-1 JPEG压缩色度和亮度量化表
此外，由于人眼对低频分量的图像比对高频分量的图像更敏感，因此表14-1中的左上角的量化步距要比右下角的量化步距小。
亮度和色度因为代表的图像的信息量不同，亮度代表了图像的低频分量，色度代表了图像的高频分量，要分别对亮度和色度进行量化，所以量化表也是不同的。
JPEG压缩色度和亮度量化表如表14-1所示。
量化会产生误差，表14-1综合了大量的图像测试的实验结果，对于大部分图像都有很好的结果。由表14-1中可以看出，高频部分对应的量化值大，目的就是将高频部分编程接近于0，以便以后处理。JPEG可以在压缩比和图像质量间作取舍，方法就是改变量化值。如果量化值放大一倍，则有更多的系数量化为0，提高了压缩比。
2.量化系数的编排
经过DCT变换后，低频分量集中在左上角，其中F(0，0)(即第一行第一列元素)代表了直流(DC)系数，即8×8子块的平均值，要对它单独编码。
由于两个相邻的8×8子块的DC系数相差很小，所以对它们采用差分编码DPCM，可以提高压缩比，也就是说，对相邻的子块DC系数的差值进行编码。8×8的其他63个元素是交流(AC)系数，采用行程编码。所以量化后的系数要重新编排，目的是为了增加连续的“0”系数的个数，就是“0”的游程长度，方法是按照Z字形的式样编排。
DCT变换后低频分量多呈圆形辐射状向高频率衰减，因此可以看成按Z字形衰减。因此，量化系数按Z字形扫描读数，这样就把一个8×8的矩阵变成一个1×64的矢量，频率较低的系数放在矢量的顶部。量化后的DCT系数的编排如图14-5所示。
图14-5 量化DCT系数的编排
量化后的DCT系数的序号如表14-2所示。
表14-2 量化后的DCT系数的序号
3.组成位数据流
JPEG编码的最后一个步骤是把各种标记代码和编码后的图像数据组成一帧一帧的数据，这样做的目的是为了便于传输、存储和译码器进行译码，这样组织的数据通常称为JPEG位数据流(JPEG Bit Stream)。
举个例子来说明上述过程。下面为8×8的亮度(Y)图像子块经过量化后的系数，如表14-3所示。
表14-3 亮度量化后的系数
可见量化后只有左上角的几个点(低频分量)不为零，这样采用行程编码就会很有效。
第一步，熵编码的中间格式表示。
先看DC系数。假设前一个8×8子块DC系数的量化值为12，则本块DC系数与它的差为3，如表14-4所示。
表14-4 DC系数表
查表得Size=2，Amplitude=3，所以DC中间格式为(2)(3)。
下面对AC 系数编码。经过Zig-Zag 扫描后，遇到的第一个非零系数为−2，其中遇到零的个数为1(即RunLength)，根据下面这张AC系数表如表14-5所示。
表14-5 AC系数表
查表得Size=2。所以RunLength=1，Size=2，Amplitude=3，所以AC中间格式为(1，2)(−2)。
其余的点类似，可以求得这个8×8子块熵编码的中间格式为：
(DC)(2)(3)，(1，2)(−2)，(0，1)(−1)，(0，1)(−1)，(0，1)(−1)，(2，1) (−1)，(EOB)(0，0)。
表14-6 AC亮度Huffman表
第二步，熵编码。
对于(2)(3)：2查表14-4得到11，3经过VLI编码为011；
对于(1，2)(−2)：(1，2)查表14-6得到11011，−2是2的反码，为01；
对于(0，1)(−1)：(0，1)查表14-6得到00，−1是1的反码，为0；……
最后，这一8×8子块亮度信息压缩后的数据流为11011，1101101，000，000，000，111000，1010。总共31比特，其压缩比是64×8/31≈16:1，因此大约每个像素用1个半比特(24/16=1.5)。
以上是图像压缩编码的整个过程。
JPEG格式是目前网络上最流行的图像格式，是可以把文件压缩到最小的格式，在Photoshop软件中以JPEG格式储存时，提供11级压缩级别，以0～10级表示。其中0级压缩比最高，图像品质最差。即使采用细节几乎无损的 10 级质量保存时，压缩比也可达5:1。以BMP格式保存时得到4.28MB图像文件，在采用JPEG格式保存时，其文件仅为178格式，压缩比达到24:1。经过多次比较，采用第8 级压缩为存储空间与图像质量兼得的最佳比例。
JPEG格式的应用非常广泛，特别是在网络和光盘读物上，都能找到它的身影。目前各类浏览器均支持JPEG这种图像格式，因为JPEG格式的文件尺寸较小，下载速度快。
14.2 基于DCT的图像压缩MATLAB仿真实现
基于DCT的JPEG图像压缩编码理论算法做应用开发仍是工程领域中的主流，因此有必要进行简要介绍。
14.2.1 数字图像文件的读写
MATLAB为用户提供了专门的函数来从图像格式的文件中读写图像数据。这种方法不像其他编程语言，需要编写复杂的代码，只需要简单地调用MATLAB提供的函数即可。
(1)图像文件的读取。MATLAB中利用函数imread 来实现图像文件的读取操作。其调用格式为：
A=imread(filename，fmt)
[X，map]=imread(filename，fmt)
[…]=imread(filename)
其中：Filename是图像文件名；fmt是图像文件格式。
A=imread(filename，fmt) 读取图像到A，如果文件是包含灰度图像，A是二维矩阵；如果文件是包含真彩色图像(RGB)，A是三维矩阵(M-by-N-by-3)。文件必须在当前目录下，或在MATLAB的路径上。如果imread不能够找到名称为filename的文件，那么它将找一名为filename. fmt的文件。
[X，map]=imread (filename，fmt) 把图像 filename 读入与它相关的图像色彩信息写入map，图像色彩信息值在范围[0，1]中自动地重新调整。
[…]=imread(filename)这种方式是试图从文件所包含的信息中得到文件的格式。
(2)图像文件的输出。MATLAB中利用函数imwrite 来实现图像文件的输出和保存操作。其调用格式格式为：
imwrite(A，filename，fmt)
imwrite(X，map，filename，fmt)
imwrite(…，filename)
imwrite(A，filename，fmt)%把图像 A 写入图像文件 filename
imwrite(X，map，ilename，fmt)%把 X和它的相关色彩信息map 写入filename
imwrite(…，filename)把图像%写入图像文件filename，并推测可能的格式用来做filename的扩展名，扩展名必须是fmt中一合法名。
(3)图像文件的显示。在MATLAB中，可以调用image函数来显示图像，其调用格式格式为：
image(C)
image(C) 把矩阵 C 转成一图像。C 可以是一M×N或M×N×3 维的矩阵，且可以包含 double、uint8或uint16 数据。image 用来显示附标图像，即显示的图像上有 x、y坐标轴的显示，可以看到图像的像素大小。但可以加上 axis off命令即可把坐标去掉。
在MATLAB的图像处理工具箱中还提供了一个应用很广泛的图像显示函数，即imshow函数，其调用格式为：
imshow(I，n)
imshow(I，[low high])
imshow(BW)
imshow(X，map)
imshow(RGB)
imshow filename
imshow(I，n)%显示一幅n个灰度级的图像I
imshow(I，[low high])%一定灰度范围内显示灰度图像，low和high参数分别为数据数组的最小值和最大值
imshow(BW)%显示二进制图像
imshow(X，map)%显示索引图像，其中X代表索引图像的数据矩阵，map为颜色映射表
imshow filename %显示一幅在当前目录下的合法文件
imshow(RGB)%显示RGB图像。RGB是一个m×n×3的数组。对于RGB中的每个像素(r，c)，imshow显示数值(r，c，1:3)所描述的颜色
例如下面的程序：
imread('zyj.bmp');
imshow(I)，title('图像读出')
14.2.2 程序流程图
基于DCT的图像压缩编码的MATLAB实现的流程图如图14-6所示。
14.2.3 DCT变换的编程实现
DCT变换是数字图像处理中重要的变换，很多重要的图像算法、图像应用都是基于DCT变换的，如JPEG图像编码方式。对于大尺寸的二维数值矩阵，倘若采用普通的DCT变换来进行，其所花费的时间将是让人难以忍受甚至无法达到实用的。而要克服这一难点，DCT变换的快速算法是比较好的选择。
就目前而言，DCT变换的快速算法无非有以下两种方式。
(1)由于FFT算法的普遍采用，直接利用FFT来实现DCT变换的快速算法相比来说就相对容易。但是此种方法也有不足：计算过程会涉及复数的运算。
由于DCT变换前后的数据都是实数，计算过程中引入复数，而一对复数的加法相当于两对实数的加法，一对复数的乘法相当于四对实数的乘法和两对实数的加法，显然是增加了运算量，也给硬件存储提出了更高的要求。
图14-6 程序流程图
(2)直接在实数域进行DCT快速变换。显然，这种方法相比于前一种而言，计算量和硬件要求都要优于前者。
在MATLAB的图像处理工具箱中，可以直接调用dct2和idct2来实现二维离散余弦变换及其反变换。
(1)dct2。
dct2函数实现图像的二维离散余弦变换，其调用格式为：
F=dct2(f)
【例14.1】将给定图像利用DCT进行变换。MATLAB代码设置如下：
f=imread('cameraman.tif');
f=im2double(f);
F=dct2(f);
subplot(121)，imshow(f，[ ]);
subplot(122)，imshow(log(1+20*abs(F))，[ ]);
运行结果如图14-7所示。
图14-7 图像显示及图像DCT变换后频谱显示
由运行结果可知，DCT变换具有能量集中的性质，数据集中在左上角。因此进行图像压缩时，离散余弦变换矩阵可以舍弃右下角的高频数据。
(2)idct2。
idct2函数实现图像的二维离散余弦逆变换，其调用格式为：
F=idct2(f)
在MATLAB图像处理工具箱中，有一个对图像进行块操作的函数blkproc，利用这个函数可以直接实现图像一系列8×8子块的DCT变换。其调用格式格式为：
B=blkproc(A，[m n]，fun，parameter1，parameter2，...)
B=blkproc(A，[m n]，[mborder nborder]，fun，...)
B=blkproc(A，'indexed'，...)
[m n]是指图像以m×n为分块单位，对图像进行处理(如8像素×8 像素)。
Fun：应用此函数分别对每个m×n分块的像素进行处理
parameter1，parameter2：要传给fun函数的参数。
mborder nborder：对每个 m×n 块上下进行 mborder个单位的扩充，左右进行 nborder个单位的扩充，扩充的像素值为0。fun函数对整个扩充后的分块进行处理。
【例14.2】利用DCT变换对灰度图实现其重构。MATLAB代码设置如下：
RGB=imread('autumn.tif');
I=rgb2gray(RGB);
J=dct2(I);
imshow(log(abs(J))，[])，
colormap(jet)，
colorbar
运行结果如图14-8所示。
J(abs(J)<10)=0;
K=idct2(J);
figure，
imshow(I)
运行结果如图14-9所示。
图14-8 色彩分布图
图14-9 原始图
figure，
imshow(K，[0 255])
运行结果如图14-10所示。
图14-10 经过反条弦变换重构后的图
14.3 基于小波压缩函数进行图像压缩
小波分析用于图像压缩的基本思想就是把图像进行多分辨率分解，分解成不同空间、不同频率的子图像，然后再对子图像进行系数编码。
14.3.1 小波变换压缩函数的应用实例
系数编码是小波变换用于压缩的核心，压缩的实质是对系数的量化压缩。图像经过小波变换后生成的小波图像的数据总量与原图像的数据量相等，即小波变换本身并不具有压缩功能。
之所以将它用于图像压缩，是因为生成的小波图像具有与原图像不同的特性，表现在图像的能量主要集中于低频部分，而水平、垂直和对角线部分的能量则较少；水平、垂直和对角线部分表征了原图像在水平、垂直和对角线部分的边缘信息，具有明显的方向特性。低频部分可以称为亮度图像，水平、垂直和对角线部分可以称为细节图像。
对所得的子图，根据人类的视觉生理和心理特点分别作不同策略的量化和编码处理。人眼对亮度图像部分的信息特别敏感，对这一部分的压缩应尽可能减少失真或者无失真。一个图像作小波分解后，可得到一系列不同分辨率的子图像，不同分辨率的子图像对应的频率是不同的。
高分辨率(高频)子图像上大部分点的数值都接近于 0，分辨率越高越明显。而对于一个图像来说，表现图像的最主要的部分是低频部分，所以最简单的压缩方法是利用小波分解去掉图像的高频部分而只保留低频部分。
小波变换中在MATLAB所使用的函数与在前面的章节中都有所介绍，这里就不再一一举例说明了。
【例14.3】使用全局阈值对给定的图像进行压缩处理。MATLAB代码设置如下：
load wmandril;
%装入待压缩图像
nbc=size(map，1);
colormap(gray(nbc));
image(wcodemat(X，nbc));
title('原图像');
axis square;
%显示原始图像
结果如图14-11所示。
[C，S]=wavedec2(X，2，'db4');
%对图像进行小波分解
thr=20;
%设置小波系数阈值
[Xcompress1，cxd，lxd，perf0，perfl2]=wdencmp('gbl'，C，S，'db4'，2，thr，'h'，2);
%使用wdencmp函数，全局阈值进行压缩，此时对所有的高频系数进行相同的阈值化处理
figure;
image(wcodemat(Xcompress1，nbc));
colormap(gray(nbc));
title(['压缩图像阈值='，num2str(thr)])
axis square
%显示压缩图像
结果如图14-12所示。
图14-11 原图
图14-12 压缩图
disp('小波系数中置0的系数个数百分比:')
perfl2
disp('压缩后图像剩余能量百分比:')
perf0
结果为：
小波系数中置0的系数个数百分比:
perfl2=
99.6127
压缩后图像剩余能量百分比:
perf0=
57.4769
【例14.4】对在水平、垂直、对角3个方向使用层相关阈值给定的图像进行压缩处理。MATLAB代码设置如下：
load wmandril;
%装入待压缩图像
subplot(1，2，1);
nbc=size(map，1);
colormap(gray(nbc));
image(wcodemat(X，nbc));
title('原始图像')
axis square
%显示原图像
thr_h=[16 19];
%水平阈值
thr_d=[20 21];
%对角阈值
thr_v=[22 23];
%垂直阈值
thr=[thr_h; thr_d; thr_v];
[Xcompress2，cxd，lxd，perf0，perfl2]=wdencmp('lvd'，X，'db3'，2，thr，'h');
%使用wdencmp函数进行压缩，在水平、垂直、对角3个方向使用层相关阈值
%保留图像小波分解的近似系数
subplot(1，2，2);
image(wcodemat(Xcompress2，nbc));
title('压缩图像，阈值=[16 19;20 21;22 23]')
axis square
%显示压缩图像
disp('小波系数中置0的系数个数百分比: ')
perfl2
disp('压缩后图像剩余能量百分比: ')
perf0
结果如图14-13所示。
图14-13 结果图
结果为：
小波系数中置0的系数个数百分比:
perfl2=
99.6099
压缩后图像剩余能量百分比:
perf0=
55.0541
【例14.5】利用小波分解去掉给定图像的高频部分而只保留低频部分。MATLAB代码设置如下：
clear
clc
load woman;%装入图像
subplot(2，2，1);
image(X);
colormap(map);
title('原始图像');
%显示原始图像
disp('压缩前图像X的大小:');
whos('X');
[c，s]=wavedec2(X，2，'bior3.7');
ca1=appcoef(c，s，'bior3.7'，1);
subplot(2，2，2);
image(ca1);
colormap(map);
title('分解后的低频系数和高频系数')
%提取小波分解结构中一层的低频系数和高频系数
ch1=detcoef2('h'，c，s，1);
%水平方向(垂直方向和斜线方向与此类似)
h1=wrcoef2('h'，c，s，'bior3.7'，1);
%重构一层水平分量(垂直方向和斜线方向与此类似)
ca1=appcoef2(c，s，'bior3.7'，1);
ca1=wcodemat(ca1，440，'mat'，0);
ca2=appcoef2(c，s，'bior3.7'，2);
ca2=wcodemat(ca2，440，'mat'，0);
subplot(2，2，3);
image(ca1);
colormap(map);
title('第一次压缩图像的大小为：');
subplot(2，2，4);
image(ca2);
colormap(map);
title('第二次压缩图像的大小为：');
结果如图14-14所示。运行结果为：
图14-14 结果图
压缩前图像X的大小:
Name　　Size　　　Bytes Class　Attributes
X　　256x256　　　524288　double
第一次压缩图像的大小为:
Name　　Size　　　Bytes Class　Attributes
ca1　　135x135　　　145800　double
第二次压缩图像的大小为:
Name　　Size　　Bytes Class　Attributes
ca2　　75x75　　　45000　double
14.3.2 基于小波包变换的图像压缩
下面我通过一个例子来说明小波包分析在图像压缩中的应用，并给出性能参数以便于同基于小波分析的压缩进行比较。
【例14.6】利用小波包分析进行图像压缩。MATLAB代码设置如下：
%读入图像
load julia
%求颜色索引表长度
nbc=size(map，1);
%得到图像的阈值，保留层数，小波树优化标准
[thr，sorh，keepapp，crit]=ddencmp('cmp'，'wp'，X)
%通过以上得到的参数对图像进行压缩
[xd，treed，perf0，perfl2]=wpdencmp(X，sorh，4，'sym4'，crit，thr*2，keepapp);
%更改索引表为pink索引表
colormap(pink(nbc));
结果如图14-15所示。
subplot(121);
image(wcodemat(X，nbc));
title('原始图像');
subplot(122);
image(wcodemat(xd，nbc));
title('全局阈值化压缩图像');
xlabel(['能量成分'，num2str(perfl2)，'%'，'零系数成分'，num2str(perf0)，'%']);
plot(treed);
结果如图14-16所示。结果为：
thr=
5.0000
sorh=
h
keepapp=
1
crit=
threshold
图14-15 结果图
图14-16 压缩图
14.4 MATLAB二维小波工具箱压缩
前面的章节曾经向读者介绍过如何利用小波工具箱的界面来进行对图像去噪的操作，本节主要介绍关于图像压缩的工具箱。
14.4.1 二维离散小波界面式压缩
本小节将向读者介绍利用小波工具箱的二维离散小波变换去除图像的噪声。
打开二维小波分析工具。首先，在窗口中键入wavemenu，弹出小波工具箱主界面。下面向读者介绍二维离散小波分析工具分析图像压缩的具体操作步骤，具体操作步骤如下。
(1)打开二维小波分析工具。首先，在窗口中键入wavemenu，弹出小波工具箱主界面，然后单击Wavelet 2-D按钮。
(2)下载图像源。选择File→Load→Signal命令，在弹出的对话框中选择woman.mat文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样图像就被下载到了此工具中，如图14-17所示。
(3)执行小波系数分解。在图14-17所示的界面中选择db4小波进行3级分解。同时单击Analyze按钮，经过短暂的计算后，工具显示分解系数图，如图14-18所示。