数据压缩实验指导书

《数据压缩》实验指导书（24学时）

实验一 图像输入与输出基本操作（2学时） 实验二 图像预测编码（3学时） 实验三 图像熵编码与压缩（3学时） 实验四 图像DCT变换编码与压缩（4学时） 实验五 图像小波变换编码与压缩（4学时） 实验六 JPEG编解码（4学时） 实验七 矢量量化图像压缩（4学时）

实验一 图像输入与输出基本操作

一、实验题目：

图像输入与输出操作

二、实验目的

学习在MATLAB环境下对图像文件的I/O操作，为读取各种格式的图像文

件和后续进行图像处理打下基础。

三、实验内容

利用MATLAB为用户提供的专门函数从图像格式的文件中读/写图像数据、显示图像，以及查询图像文件的信息。 四、预备知识

熟悉MATLAB开发环境。

五、实验原理

（1）图像文件的读取

利用imread函数可以完成图像文件的读取操作。常用语法格式为： I=imread(‘filename’,‘fmt’)或I=imread(‘filename.fmt’)；

其作用是将文件名用字符串filename表示的、扩展名用字符串fmt（表示图像文件格式）表示的图像文件中的数据读到矩阵I中。当filename中不包含任何路径信息时，imread会从当前工作目录中寻找并读取文件。要想读取指定路径中的图像，最简单的方法就是在filename中输入完整的或相对的地址。MATLAB支持

多种图像文件格式的读、写和显示。因此参数fmt常用的可能值有：

‘bmp’ Windows位图格式 ‘jpg’or‘jpeg’ 联合图像专家组格式 ‘tif’or‘tiff’ 标志图像文件格式 ‘gif’ 图形交换格式 ‘pcx’ Windows画刷格式 ‘png’ 可移动网络图形格式 ‘xwd’ X Window Dump格式 例如，命令行 >>I=imread(‘lena.jpg’);

将JPEG图像lena读入图像矩阵I中。 （2） 图像文件的写入（保存）

利用imwrite完成图像的输出和保存操作，也完全支持也完全支持上述各种 图像文件的格式。其语法格式为：

imwrite(I,‘filename’,‘fmt’)或imwrite(I,‘filename.fmt’)； 其中的I、filename和fmt的意义同上所述。

注意事项：当利用imwrite函数保存图像时，MATLAB默认的保存方式是将其简化为uint8的数据类型。与读取文件类型类似，MATLAB在文件保存时还支持16位的PNG和TIFF图像。所以，当用户保存这类文件时，MATLAB就将其存储在uint16中。 （3）图像文件的显示

图像的现实过程是将数字图像从一组离散数据还原为一幅可见图像的过程。 MATLAB的的图像处理工具箱提供了多种图像显示技术。例如imshow可以直接从文件显示多种图像；image函数可以将矩阵作为图像 ；colorbar函数可以用来显示颜色条；montage函数可以动态显示图像序列。这里仅对常用的显示函数进行介绍。

①图像的显示

imshow函数是最常用的显示各种图像的函数，其调用格式如下：

imshow(I,N)；

imshow(I,N)用于显示灰度图像，其中I为灰度图像的数据矩阵，N为灰度级

数目，默认值为256。

例如下面的语句用于显示一幅灰度图像： >> I=imread(‘lena.jpg’)； >> imshow(I);

如果不希望在显示图像之前装载图像，那么可以使用以下格式直接进行图像

文件的显示：

imshow filename

其中，filename为要显示的图像文件的文件名。 实例1-1 显示一幅在当前目录下的.bmp格式的图像： >>imshow lena.jpg 显示结果如图1.1所示。 图1.1.1 显示一幅图像文件中的图像

注意事项：该文件名必须带有合法的扩展名（指明文件格式），且该图像文件必须保存在当前目录下，或在MATLAB默认的目录下。

②添加色带

colorbar函数可以给一个坐标轴对象添加一条色带。如果该坐标轴对象包含一个图像对象，则添加的色带将指示出该图像中不同颜色的数据值。这对于了解被显示图像的灰度级特别有用。其调用格式为：

colorbar

实例1-2

>> I=imread(‘lena.jpg’)； >> imshow(I); >> colorbar;

图1.1.2 显示图像并加入颜色条

从上图可知，该图像是数据类型为uint8的灰度图像，其灰度级范围从0 －255。

③显示多幅图像

显示多幅图像最简单的方法就是在不同的图形窗口中显示它们。imshow总

是在当前窗口中显示一幅图像，如果用户想连续显示两幅图像，那么第二幅图像就会替代第一幅图像。为了避免图像在当前窗口中的覆盖现象，在调用imshow函数显示下一幅图像之前可以使用figure命令来创建一个新的窗口。例如：

imshow(I1); figure, imshow(I2); figure, imshow(I3);

有时为了便于在多幅图像之间进行比较，需要将这些图像显示在一个图形

窗口中。达到这一目的有两种方法：一种方法是联合使用imshow和subplot函数，但此方法在一个图形窗口只能有一个调色板；另一种方法是联合使用subimage和subplot函数，此方法可在一个图形窗口内使用多个调色板。

subplot函数将一个图形窗口划分为多个显示区域，其调用格式如下：

subplot(m,n,p)

subplot函数将图形窗口划分为m（行）×n（列）个显示区域，并选择第p个区域作为当前绘图区。

例1-3 用两排显示四幅图像，可以使用以下语句： >> I1=imread(‘lena.bmp’)；

>> I2=imread(‘gs256.bmp’)； >> I3=imread(‘lena.bmp’)； >> I4=imread(‘gs256.bmp’)； >> subplot（2,2,1）, imshow(I1); >> subplot（2,2,2）, imshow(I2); >> subplot（2,2,3）, imshow(I3); >> subplot（2,2,4）, imshow(I4);

图1.1.3 在一个图形窗口中显示多幅图像

（4） 图像文件信息的查询

imfinfo函数用于查询图像文件的有关信息，详细地显示出图像文件的各种属性。其语法格式为：

info=imfinfo(‘filename’,‘fmt’)或info=imfinfo(‘filename.fmt’) 或imfinfo filename.fmt

imfinfo函数获取的图像文件信息依赖于文件类型的不同而不同，但至少应包 含以下内容：

文件名。如果该文件不在当前目录下，还包含该文件的完整路径。 文件格式。 文件格式的版本号。 文件最后一次修改的时间。

文件的大小。以字节为单位。 图像的宽度。 图像的高度。

每个像素所用的比特数。也叫像素深度。

图像类型。即该图像是真彩色图像、索引图像还是灰度图像。 例如，命令行

>>imfinfo bubbles25.jpg

会输出如下信息（注意，在这种情况下，有些域不包含信息）：

Filename: ‘bubbles25.jpg’

FileModDate: ‘04-Jan-2003 12:31:26’

FileSize: 13849 Format: ‘jpg’

FormatVersion: ‘’

Width: 714 Height: 682 BitDepth: 8

ColorType: ‘grayscale’ FormatSignature: ‘’

Comment: {}

六、实验步骤：

（1）利用imread函数完成对图像文件的读取操作。 （2）利用imwrite函数完成图像的写入（保存）操作。 （3）利用imshow函数显示图像。

（4）利用imfinfo函数查询图像文件的有关信息。

七、思考题目：

（1）若分别读入空间分辨率为M×N的灰度图像和彩色图像，则得到的矩阵

的维数是否相同？如果不同，则分别是几维的矩阵？ （2） 读入一幅索引图像。 （3） 动态显示一个图像序。

实验二 图像预测编码

一、实验题目:

图像预测编码：

二、实验目的:

实现图像预测编码和解码.

三、实验内容:

给定一幅图片，对其进行预测编码，获得预测图像，绝对残差图像, 再利用预测图像和残差图像进行原图像重建并计算原图像和重建图像误差.

四、预备知识:

预测方法，图像处理概论。

五、实验原理:

根据图像中任意像素与周围邻域像素之间存在紧密关系,利用周围邻域四个像素来进行该点像素值预测,然后传输图像像素值与其预测值的差值信号，使传输的码率降低，达到压缩的目的。

六、实验步骤:

(1) 选取一幅预测编码图片;

(2) 读取图片内容像素值并存储于矩阵; (3) 对图像像素进行预测编码; (4) 输出预测图像和残差图像; (5) 根据预测图像和残差图像重建图像; (6) 计算原预测编码图像和重建图像误差.

七、思考题目:

如何采用预测编码算法实现彩色图像编码和解码.

八、实验程序代码:

预测编码程序1:

%编码程序:

i1=imread('lena.jpg');

if isrgb(i1) %判断i1是否为reb图像 i1=rgb2gray(i1); %如果i1为rgb图像，那么将其转换为灰度图像 end

i1=imcrop(i1,[1 1 256 256]); %细节放大[a,b,c,d] a,b为x，y的起始

%位置，c,d为x,y的宽度

i=double(i1);

[m,n]=size(i); %m为矩阵i的行大小、n为矩阵i的列大小 p=zeros(m,n); %将p置为m*n的零矩阵 y=zeros(m,n); %将y置为m*n的零矩阵

y(1:m,1)=i(1:m,1); %将矩阵i的第1列复制给y矩阵的第1列 p(1:m,1)=i(1:m,1); %将矩阵i的第1列复制给p矩阵的第1列 y(1,1:n)=i(1,1:n); %将矩阵i的第1行复制给y矩阵的第1行 p(1,1:n)=i(1,1:n); %将矩阵i的第1行复制给p矩阵的第1行 y(1:m,n)=i(1:m,n); %将矩阵i的第n列复制给y矩阵的第n列 p(1:m,n)=i(1:m,n); %将矩阵i的第n列复制给p矩阵的第n列 p(m,1:n)=i(m,1:n); %将矩阵i的第m行复制给y矩阵的第m行 y(m,1:n)=i(m,1:n); %将矩阵i的第m行复制给p矩阵的第m行 for k=2:m-1

for l=2:n-1

y(k,l)=(i(k,l-1)/2+i(k-1,l)/4+i(k-1,l-1)/8+i(k-1,l+1)/8); p(k,l)=round(i(k,l)-y(k,l)); %round为四舍五入取整

%p为原图像与编码后图像的差值

end end

p=round(p); %将p内所有的值取整 subplot(3,2,1); imshow(i1); title('原灰度图像'); subplot(3,2,2);

imshow(y,[0 256]); %显示y的灰度范围从0到256 title('利用三个相邻块线性预测后的图像'); subplot(3,2,3); imshow(abs(p),[0 1]); title('编码的绝对残差图像');

%解码程序 j=zeros(m,n); j(1:m,1)=y(1:m,1); j(1,1:n)=y(1,1:n); j(1:m,n)=y(1:m,n); j(m,1:n)=y(m,1:n); for k=2:m-1

for l=2:n-1

j(k,l)=p(k,l)+y(k,l); end end for r=1:m

for t=1:n

d(r,t)=round(i1(r,t)-j(r,t)); end end

subplot(3,2,4); imshow(abs(p),[0 1]); title('解码用的残差图像'); subplot(3,2,5); imshow(j,[0 256]);

title('使用残差和线性预测重建后的图像'); subplot(3,2,6); imshow(abs(d),[0 1]);

title('解码重建后图像的误差');

九、实验结果:

图1.10.1 Lena图像预测编码实验结果

预测编码程序2:

x=imread('e:\\imagebase\\cameraman.jpg'); figure(1) subplot(2,3,1); imshow(x); title('原始图像'); subplot(2,3,2); imhist(x);

title('原始图像直方图'); subplot(2,3,3); x=double(x); x1=yucebianma(x); imshow(mat2gray(x1));

title('预测误差图像'); subplot(2,3,4); imhist(mat2gray(x1)); title('预测误差直方图'); x2=yucejiema(x1); subplot(2,3,5); imshow(mat2gray(x2)); title('解码图像'); e=double(x)-double(x2); [m,n]=size(e);

erms=sqrt(sum(e(:).^2)/(m*n)); %平方和开方

%预测编码函数；

%一维无损预测编码压缩图像x，f为预测系数，如果f默认，则f=1，即为前值预测

function y=yucebianma(x,f) error(nargchk(1,2,nargin)) if nargin<2 f=1; end

x=double(x); [m,n]=size(x); p=zeros(m,n); xs=x;

zc=zeros(m,1); if length(f)>1 for j=1:length(f)

xs=[zc xs(:,1:end-1)]; p=p+f(j)*xs; end

y=x-round(p); else

xs=[zc xs(:,1:end-1)]; p=xs; y=x-round(p); end

% yucejiema是解码程序,与编码程序用的是同一个预测器 function x=yucejiema(y,f) error(nargchk(1,2,nargin)); if nargin<2 f=1; end

if length(f)>1 f=f(end:-1,1); [m,n]=size(y); order=length(f); x0=zeros(m,n+order); for j=1:n jj=j+order; for i=1:m tep=0.0; for k=order:-1:1

tep=tep+f(k)*x0(i,jj-order-1+k); end

x0(i,jj)=y(i,j)+tep; end end

x=x0(:,order+1:end); else

[m,n]=size(y); x0=zeros(m,n+1); for j=1:n jj=j+1; for i=1:m

x0(i,jj)=y(i,j)+x0(i,jj-1); end end

x=x0(:,2:end); end

图1.10.2 摄影师图像预测编码实验结果

实验三 图像熵编码与压缩

一、实验题目：

图像熵编码与压缩

二、实验目的：

学习和理解建立在图像统计特征基础上的熵编码压缩方法。

三、实验内容：

（1）编程实现二值文本图像的行程编码。

（2）编程实现灰度图像的霍夫曼编码，并计算图像熵、平均码字长度及

编码效率。

四、预备知识：

（1）熟悉行程编码原理。 （2）熟悉霍夫曼编码原理。

（3）熟悉在MATLAB环境下对图像文件的I/O操作。

五、实验原理：

（1）行程编码

行程编码是将一行中颜色相同的相邻像素用一个计数值和该颜色值来代替。 比如，aaabbcccccdddeeee可以表示为3a2b5c3d4e。如果一幅图像由很多块颜色相同的大面积区域组成，则采用行程编码可大大提高压缩效率，尤其适用于二值图像。但当图像中每两个相邻像素的颜色都不相同时，采用这种方法不但不能实现数据压缩，反而使数据量增加一倍。因此，对复杂的图像都不能单纯地采用行程编码。

（2） 霍夫曼编码

霍夫曼编码是一种代码长度不均匀的编码方法。它的基本原理是按信源符号出现的概率大小进行排序，出现概率大的分配短码，反之则分配长码。

霍夫曼编码基本步骤如下：

步骤1：统计图像每个灰度级（信息符号）出现的概率，并按概率从大到小进行排序。

步骤2：选出概率最小的两个值进行组合相加，形成的新概率值和其他概率 值形成一个新的概率集合。

步骤3：重复步骤2，反复利用合并和排序的方法，直到只有两个概率为止。 步骤4： 分配码字，对最后两个概率一个赋予“0”码字，一个赋予“1”码字。

如此反向进行到开始的概率排列，这样就得到了各个符号的霍夫曼编码。

六、实验步骤：

（1）编程实现二值文本图像的行程编码。

（2）编程实现连续灰度图像的霍夫曼编码，并计算图像熵、平均码字长度及编码效率。

七、思考题目：

将行程编码与霍夫曼编码结合，能否提高压缩效果？试验证之。

八、实验程序代码：

（1）二值文本图像的行程编码程序： clear,close all t=imread('text.tif'); ts=logical(t);

codetable=zeros(1,20000); [m,n]=size(ts); nn=n+1; icodecount=1; for i=1:m

p1=ts(i,1);%第i行，第1个像素 ipcount=1;%同一灰度值连续出现的次数 for j=2:n p2=ts(i,j);

if ((p1==p2)&(ipcount

codetable(icodecount)=ipcount; codetable(icodecount+1)=p1; icodecount=icodecount+2; p1=p2; ipcount=1; end;

end;

codetable(icodecount)=ipcount; codetable(icodecount+1)=p1; icodecount=icodecount+2; codetable(icodecount)=nn;%行结束符号 icodecount=icodecount+1; end;

codetable(icodecount)=65535;%码表结束符号 （2）连续灰度图像的霍夫曼编码程序代码 function s=reduce(p) s=cell(length(p),1); for i=1:length(p) s{i}=i; end;

n=size(s,1); while n>2

[p,i]=sort(p); p(2)=p(1)+p(2); p(1)=[]; s=s(i);

s{2}={s{1},s{2}}; s(1)=[];

n=size(s,1); end;

function makecode(sc,codeword) global CODE if isa(sc,'cell')

makecode(sc{1},[codeword 0]); makecode(sc{2},[codeword 1]); else

CODE{sc}=char('0'+codeword); end;

function CODE=huffman(p) global CODE

CODE=cell(length(p),1); if length(p)>1 p=p/sum(p); s=reduce(p); makecode(s,[]);

else

CODE={'1'}; end;

function y=mat2huff(x) y.size=uint32(size(x)); x=round(double(x)); xmin=min(x(:)); xmax=max(x(:));

pmin=double(int16(xmin)); pmin=uint16(pmin+32768); y.min=pmin;

x=x(:)';

h=histc(x,xmin:xmax); maxh=max(h); if maxh>65535

h=65535*h/maxh; end;

h=uint16(h); y.hist=h;

map=huffman(double(h)); hx=map(x(:)-xmin+1); hx=char(hx)'; hx=hx(:)'; hx(hx==' ')=[];

ysize=ceil(length(hx)/16); hx16=repmat('0',1,ysize*16); hx16(1:length(hx))=hx;

hx16=reshape(hx16,16,ysize); hx16=hx16'-'0';

twos=pow2(15:-1:0);

%y.code=uint16(sum(t3,2))'; t1=ones(ysize,1); t2=twos(t1,:); t3=hx16.*t2; t4=sum(t3,2);

y.code=uint16(t4)';

实验二 图像DCT变换编码与压缩

一、实验题目：

图像DCT变换编码与压缩

二、实验目的：

（1）掌握离散余弦变换DCT的实现方法，了解DCT的幅度分布特性，从而加深对DCT变换的认识。

（2）掌握图像DCT变换编码的实现方法，从而加深对变换编码压缩图像原理的理解。

三、实验内容：

编程实现图像DCT变换编码。

四、预备知识：

（1）熟悉DCT原理。 （2）熟悉变换编码原理。

（3）熟悉在MATLAB环境下对图像文件的I/O操作。

五、实验原理：

变换编码是在变换域进行图像压缩的一种技术。图2.2.1显示了一个典型的

变换编码系统。

输入图像N×N构建n×n子图像正变换量化器符号编码器压缩图像

图2.2.1 变换编码系统

在变换编码系统中，如果正变换采用DCT变换就称为DCT变换编码系统。 DCT用于把一幅图像映射为一组变换系数，然后对系数进行量化和编码。对于大多数的正常图像来说，多数系数具有较小的数值且可以被粗略地量化（或者完全抛弃），而产生的图像失真较小。

DCT是仅次于K-L变换的次最佳正交变换，且以获得广泛应用，并成为许多图像编码国际标准的核心。离散余弦变换的变换核为余弦函数，计算速度快，有利于图像压缩和其他处理。

对于N×N的数字图像，二维DCT变换的正反变换可表示为：

F(u,v)?c(u)c(v)??f(x,y)cosx?0y?0N?1N?1(2x?1)u?(2y?1)v?cos2N2N(2x?1)u?(2y?1)v?cos2N2Nf(x,y)?2N??c(u)c(v)F(u,v)cosu?0v?0N?1N?1（2.2.1）

其中，

?1/2u?0或v?0?c(u)?c(v)??

??1u,v?1,2,...,N?1MATLAB图像处理工具箱实现离散余弦变换有两种方法：

（1）使用函数dct2，该函数用一个基于FFT的算法来提高当输入较大的方阵时的计算速度。

（2）使用由dctmtx函数返回的DCT变换矩阵，这种方法较适合于较小的输入方阵（例如8×8或16×16）。

①函数：dct2

实现图像的二维离散余弦变换。调用格式为： B = dct2(A) B = dct2(A,[M N]) B = dct2(A,M,N)

式中A表示要变换的图像，M和N是可选参数，表示填充后的图像矩阵大小，B表示变换后得到的图像矩阵。

②函数：dctmtx

除了用dct2函数实现二维离散余弦变换，还可用 dctmtx函数来计算变换矩阵，调用格式为：

D = dctmtx(N)

式中D是返回N×N的DCT变换矩阵，如果矩阵A是N×N方阵，则A的DCT变换可用D×A×D’来计算。这在有时比dct2计算快，特别是对于A很大的情况。

③函数：idct2

实现图像的二维离散余弦反变换。调用格式为： B = idct2(A) B = idct2(A,[M N]) B = idct2(A,M,N)

式中参数同dct2。

此外，为了实现8×8子块的DCT图像变换还要用到MATLAB中的blkproc函数。将这个函数和函数dctmtx一起用于块处理可以大大简化运算。调用函数blkproc的格式为：

B=blkpro(A,[M,N],FUN,P1,P2,…)

其中，A表示原图像，[M,N]指定了大小为M×N的滑动邻域，FUN是对M×N的矩阵进行计算的函数，Pi是传递给FUN的附加参数。该函数自动实现图像块处理的整个过程。Blkproc把A分成M×N个块，对每个块调用参数为P1，P2，…的函数FUN，并重新将结果组合到输出图像B。

blkproc函数实现n×n矩阵的DCT变换和反变换。编程中可写成： Y=blkproc(F，[8 8]，’P1*x*P2’，H，H’) 同样的道理，blkproc函数还用于量化和反量化。

显示误差直方图可能用到的MATLAB函数有： Max %找图像差最大值 [ ]=hist %用于生成直方图数据 Bar %显示图像差值直方图

以上函数用MATLAB的help查看具体使用方法。

图2.2.2显示了采用JPEG标准化矩阵进行DCT变换编码的结果。

图2.2.2 DCT变换编码

六、实验步骤：

DCT变换编码流程如下:

步骤1：设置JPEG标准化数组; 步骤2：求8×8快的DCT变换矩阵; 步骤3: 计算8×8快的DCT变换； 步骤4：对DCT系数量化和反量化； 步骤5：求反量化系数的逆DCT变换；

步骤6：重新显示重建图像、误差图像和误差图像的直方图。 量化时可采用JPEG标准推荐的归一化数组，如表2.2.1所示。

表2.2.1 JPEG标准化数组

16 12 14 14 18 24 49 72 11 12 13 17 22 35 64 92 10 14 16 22 37 55 78 95 16 19 24 29 56 64 87 98 24 26 40 51 68 81 103 112 40 58 57 87 109 104 121 100 51 60 69 80 103 113 120 103 61 55 56 62 77 92 101 99 七、思考题目： （1）观察图像8×8子块的DCT系数的分布，并分析其特点。

（2）将量化步长分别增大为初始值的2倍、4倍、8倍后再进行DCT变换编码，显示不同量化步长条件下的重建图像、误差图像以及误差图像的直方图。分析重建图像质量和量化步长的关系。

八、实验程序代码：

function y=dctcoder(x,quality) m=[16 11 10 16 24 40 51 61 12 12 14 19 26 58 60 55 14 13 16 24 40 57 69 56 14 17 22 29 51 87 80 62 18 22 37 56 68 109 103 77 24 35 55 64 81 104 113 92 49 64 78 87 103 121 120 101 72 92 95 98 112 100 103 99]*quality;

[xm,xn]=size(x); xx=x; x=double(x); t=dctmtx(8);

y1=blkproc(x,[8 8],'P1*x*P2',t,t');üT y2=blkproc(y1,[8 8],'round(x./P1)',m);%Q y3=blkproc(y2,[8 8],'x.*P1',m);%IQ y=blkproc(y3,[8 8],'P1*x*P2',t',t);%IDCT subplot(2,2,1);imshow(xx);title('原始图像');

subplot(2,2,2);imshow(mat2gray(y));title('重建图像');%reconstruted image d=x-y;%original-reconstruted

subplot(2,2,3);imshow(mat2gray(d));title('误差图像'); [h,k]=hist(d(:),512);

subplot(2,2,4);bar(k,h,'k');title('误差图像直方图');

实验五 图像小波变换编码与压缩

一、实验题目：

图像小波变换编码与压缩

二、实验目的：

掌握图像压缩编码的原理，熟悉小波变换的基本知识，了解嵌入零数小波编

码原理，并应用MATLAB编程实现嵌入零数小波编、解码程序。

三、实验内容：

（1）对图像进行小波分解求取小波系数； （2）编程实现嵌入零数小波的编码程序； （3）编程实现嵌入零数小波的解码程序。

四、预备知识：

（1）熟悉图像读写和显示； （2）理解图像压缩编码的原理； （3）理解小波变换的原理；

（4）理解嵌入零数小波编码原理。

五、实验原理：

随着多媒体和网络技术的迅猛发展,人们对基于网络的图像传输、浏览和检

索等应用提出了愈来愈广泛的需求,研究具有高效的压缩能力、支持用户个性化浏览并适合网络渐进传输的图像编码方法已成为目前静态图像编码领域的研究热点。

在大量的图像编码方法中,Shapiro提出的嵌入式零树小波编码算法

(Embedded Zerotree Wavelet encoder ，简称EZW编码)是公认的效率最高的图像渐进式编码（progressive encoding）方法之一, 这种算法得到的比特流中的比特按其重要性排序。使用这种算法，编码者能够在任一点结束编码，允许精确到任一个目标比特率或目标失真率。

（1）嵌入零数小波编码算法原理 ①内嵌编码(Embedded Coding)

内嵌编码就是编码器将待编码的比特流按重要性的不同进行排序，根据目标

码率或失真度大小要求随时结束编码；同样，对于给定码流解码器也能够随时结束解码，并可以得到相应码流截断处的目标码率的恢复图像。内嵌编码的次序是从最重要的位（最高位）到最不重要的位（最低位）逐个发送，直到达到所需码率时停止。

内嵌编码的输出信息主要包括两部分：排序信息和重要像素的位信息。位信

息是编码必不可少的有效信息；排序信息是辅助信息，按其重要性从左到右排列，反映了重要像素在原图上的空间位置，用于恢复原始的数据结构。

一幅图像经过三级小波分解后形成了10个子带，如图7-1所示。小波系数的

分布特点是越往低频子带系数值越大，包含的图像信息越多，对视觉比较重要，如图7-1中的LL3子带。越往高频子带系数值越小，包含的图像信息越少，对视觉来说不太重要。这样对相同数值的系数选择先传较低频的系数的重要比特，后传较高频系数的重要比特。正是由于小波系数具有这些特点，它非常适合于嵌入式图像的编码算法。在JPEG2000标准中以小波变换作为图像编码的变换方法。

图2.7.1 小波系数间的父子关系

②几个基本概念

一幅原始图像经过变换后得到的是空间分布的小波系数，引入字符来表示每一个小波系数状态，这样更有利于将来的量化和熵编码。

零数根(ZTR)：在本次量化中，该小波系数的幅值小于给定的阈值T，并且它的子孙后代中任何一个的幅值都比T小。与此同时，还要求该系数的父节点系数的幅值大于T；

负大系数(NEG)：在本次量化中，该小波系数的幅值大于给定的阈值T。但是，它的最高位为1，即它实际的值是一个负数；

正大系数(POS)：在本次量化中，该小波系数的幅值大于给定的阈值T。并且，它的最高位为0，即它实际的值是一个正数；

孤独零(IZ)：在本次量化中，该小波系数的幅值小于给定阈值T。但是，它的子孙后代中至少存在一个系数的幅值大于T，即该小波系数的子孙后代中存在有大系数（无论正负）。与此同时，该系数的父节点的幅值可以为任意值。

零数，在某一量化级中它的任何一个节点都不是大系数，其中，最上层的那个系数就是零数根。

EZW算法利用小波系数的特点较好的实现了图像编码的嵌入功能，主要包括：零数预测、用零数结构编码重要图和逐次逼近量化。

③零数预测

一幅经过小波变换的图像按其频带从低到高形成一个树状结构，树根是最低频子带的节点，它有3个孩子分别位于3个次低频子带的相应位置，其余子带（最高频子带除外）的节点都有4个孩子位于高一级子带的相应位置（由于高频子带

分辨率增加，所以一个低频子带节点对应有四个高频子带节点，即相邻的2×2矩阵）。这样图2.7.1所示的三级小波分解就形成了深度为4的树。

定义一个零树的数据结构：一个小波系数x，对于一个给定的门限T，如果|x|

④用零树结构编码重要图

重要图包括3种要素：重要系数、孤立零和零树根。其中，对于一个给定的阈值T，如果系数x本身和它的所有子孙都小于T，则该点就称为零树根；如果系数本身小于T，但其子孙至少有一个大于或等于T，则该点就称为孤立零点。在编码时分别用三种符号与之对应。当编码到最高分辨率层的系数时，由于它们没有子孙，零树根不再存在，只需其余两种符号即可。为了有利于内嵌编码，将重要系数的符号与重要图一起编码，这样就要使用四种符号：零树根、孤立零、正重要系数、负重要系数。

⑤逐次逼近量化（SAQ: Successive-Approximation Quantization）

内嵌编码的核心在于采用了逐次逼近的量化方法（SAQ）。SAQ按顺序使用了一系列阈值T0,T1,…,TN-1来判断重要性，其中，Ti=Ti-1/2，初始阈值T0按如下条件选择，|Xj|<2T0，其中Xj表示所有变换系数。

在编（译）码过程中，始终保持着两个分离的列表：主表和辅表。主表对应于编码中的不重要的集合或系数，其输出信息起到了恢复各重要值的空间位置结构的作用，而辅表是编码的有效信息，输出为各重要系数的二进制值。

编码分为主、辅两个过程：在主过程中，设定阈值为Ti，按上述原理对主表进行扫描编码，若是重要系数，则将其幅值加入辅表中，然后将该系数在数组中置为零，这样当阈值减小时，该系数不会影响新零树的出现；在辅过程中，对辅表中的重要系数进行细化，细化过程类似于比特平面编码。扫描顺序对最终的压缩结果会有些许影响，目前提出的扫描顺序有好几种，比较常见的有如图2.7.3所示的两种。

图2.7.3 小波系数的扫描顺序

对阈值Ti来说，重要系数的所在区间为[Ti,2Ti]，若辅表中的重要系数位于

[Ti,3Ti/2]，则用符号“0”表示，重构值为1.25 Ti-1；否则用符号“1”表示，重构值为1.75 Ti-1。输出的符号“0”、“1”由一个辅扫描表记录。

编码在两个过程中交替进行，在每个主过程前将阈值减半。译码时系数的重构值可以位于不确定区间的任意处，如果采用MMSE准则，则重构值应位于不确定区间的质心处。

SAQ是和位平面编码直接关联的，它是EZW快速算法的基础。 （2）算法分析

在图像的低比特率编码中，用来表示非零系数所在位置的开销远远大于用来表示非零系数值的开销。零树结构正是一种描述图像经过小波变换后非零数值位置的有效方法。

EZW的编码思想是不断扫描变换后的图像，生成多棵零树来对图像编码：一棵零树的形成需要对图像进行两次扫描。在生成第一棵零树时，首先找出变换后图像的最大绝对值系数，用它对应的T0作初始阈值，对图像进行第一次扫描。将图像中绝对值小于阈值的系数都看作零，然后按前面的符号定义形成零树。在第二次扫描中，对那些绝对值大于阈值的节点（POS和NEG），按其绝对值是否大于阈值的1.5倍附加一个比特1或0来描述其精度。

这样做的目的是减小非零节点系数值的变换范围，使其适应下一次阈值减半后的比特附加。而后将阈值减半，再经两次扫描生成第二棵零树，在第一次扫描

生成零树时，以前已经大于阈值的节点不再考虑，而第二次扫描附加比特时则要考虑以前数值较大的节点以保证精度。如此重复下去，不断生成零树，直到达到需要的编码精度为止。

EZW算法也存在一些问题：

①由于编码时它形成多棵零树，因而要多次扫描图像，效率低。且每一棵树必须在前一棵树形成之后才能形成，难以用并行算法优化；

②对所有频域进行等同重要度的编码，不能充分利用小波变换的特点。改进办法之一是把最低频子图与其它子图分开处理，对其进行单独得无失真编码；

③在一棵零树中包含的元素越多，越有利于数据压缩。在EZW算法中存在这样的树间冗余，在后续算法如SPIHT中则进一步利用了这种树间冗余；

④通过对小波系数的分析发现，在同一子带中相邻元素间有一定的相关性，尤其在高频子带存在大量的低值元素，所以可以通过子带中的集合把大量的这种低值元素组织到一起，达到数据压缩的目的。EZW算法并没有充分利用这种相关性。

以Lena.bmp的局部图像（128×128）为例，图2.7.4显示了该图经过三层小波变换后的小波系数分布。

图2.7.3 一幅图像三层小波变换后的小波系数分布

该图经过EZW算法编码后每级编码主表及辅表得到的部分数据流显示如下 （此处扫描顺序为RASTER scan）：

主表：Znnnznznzzzzzzznnznnnzznzzzzzzznnznnzzn

zzzzzzzznnzzzznnnzzzzzznnnzzzzzzzzzppp

辅表：0000000011001101011100011000011100000100

000000100100101101000000000

(3)嵌入零数小波解码算法

解码过程其实与主扫描过程非常类似，按照扫描次序扫描解码矩阵。 ①初始化

获取扫描次序列表，初始化符号矩阵，创建以下几个空表：重要系数重构列表，量化器编号列表，上一次解码用到的辅扫描表。

获取本级解码需要的主扫描表和辅扫描表。创建解码矩阵，量化符号列表编号，主扫描表扫描编号。

②构造逆量化器

逆量化器包括两个部分，一是量化值部分，与编码程序中的量化器函数代码相同；二是量化器编号列表的构造，这个是难点和重点，把它构造出来，才能解码出重要系数的数值。

③提高上一级解码的重要系数的重构精度

根据逆量化器生成的重构值矩阵和逆量化器编号列表t来更新重要系数列表的值，并将更新数据存入解码矩阵中。

根据不同阈值的渐进式解码后的重构图像如图2.7.4所示。

图2.7.4 不同阈值下解码后的重构图像

六、实验步骤:

（1）打开一幅灰度图像，对该图像进行三层小波变换，并得到小波系数分

布；

（2）按照EZW算法的原理实现图像的编码程序； （3）编程解码程序，实现相应的图像渐进式解码； （4）实验完毕后，提交一份实验报告。

七、思考题目:

（1）在JPEG2000标准中小波变换为什么取代了DCT变换？

（2）不同的扫描顺序对编、解码效果会有什么影响？ （3）EZW算法的编码原理是什么？ （4）EZW算法的译码过程如何？ （5）EZW算法有何缺点？该如何改进？

八、实验程序代码：

（1）生成查找表

function y=lookup_128

% Prepare table of indices of descendants y=extract(reshape([1025:2048],32,32)'); y=[y ;extract(reshape([2049:3072],32,32)')]; y=[y ; extract(reshape([3073:4096],32,32)')]; y=[y ;extract(reshape([4097:8192],64,64)')]; y=[y ; extract(reshape([8193:12288],64,64)')]; y=[y ; extract(reshape([12289:16384],64,64)')]; y=[zeros(512,4);y]; indices=y; save indices

（2）EZW编码程序

clear,clc load indices format long

% Prepare for RASTER scan [Sha93] xm=mapping_128; xm=xm(:); load lena.mat

mat=anal2d(x(50:50+127,50:50+127),3,3); x=xm; seqt(x)=mat;

T=2^round(log2(max(max(abs(mat))))); T=T/2

mat=mat.*[abs(mat)>=T==1];

seq(x)=mat; ztr=[]; sig_coeff=[]; num_of_passes=1; refine=[]; kkk=0; beta=[]; sym=[];

% For the subband LL3

while num_of_passes<=9 % Number of dominant passes kkk=kkk+1; % This is only for level 3 for ii=513:1024

if(isempty(sig_coeff) | sum([sig_coeff==ii])==0) % Added new count=0; for kk=1:4

if(seq(indices(ii,kk))==0),count=count+1; end for mm=1:4

if seq(indices(indices(ii,kk),mm))==0,count=count+1;

end

end end

if count==20 & seq(ii)==0, ztr=[ztr ii]; for kk=1:4

seq(indices(ii,kk))=inf; for mm=1:4

seq(indices(indices(ii,kk),mm))=inf; end end end end % Added new end

% Start of level 2

for ii=1025:4096

if(isempty(sig_coeff) | sum([sig_coeff==ii])==0) % Added new count=0;

if(seq(ii)~=inf & seq(ii)~=0), for kk=1:4

if(seq(indices(ii,kk))==0),count=count+1;end end

elseif seq(ii)==0 & seq(ii)~=inf, for kk=1:4

if(seq(indices(ii,kk))==0),count=count+1;end end

if count==4 & seq(ii)==0, ztr=[ztr ii]; for kk=1:4

seq(indices(ii,kk))=inf; end end end end end

%--------------------------------------------------------------------------------------- % Start of encoding for ii=1:128*128

if(isempty(sig_coeff) | sum([sig_coeff==ii])==0), if (seq(ii)>0 & seq(ii)~=inf),sym=[sym 'p']; elseif (seq(ii)<0 & seq(ii)~=inf),sym=[sym 'n']; elseif seq(ii)==0

if(ii<4097 & ~isempty(find(ztr==ii) & seq(ii)~=inf)) sym=[sym 'r'];

elseif(ii<4097 & seq(ii) ~= inf & isempty(find(ztr==ii))) sym=[sym 'z']; % This Was i earlier elseif(ii>4096 & seq(ii)~=inf) sym=[sym 'z']; end

end end end

% Encoding ends for one dominant pass

zero_tree_roots_positions=ztr; % Well just for debugging % Prepare for further passes

sig_coeff=[sig_coeff find(seq~=0 & seq~=inf)] ; % Index of where significant coeff occur

% Perform the Subordinate / Refinement pass

ref=dec2bin(abs(round(seqt(sig_coeff))),length(dec2bin(round(max(max(abs(seqt))))))); % Try not calculating the calculated % Change this,this takes long time

ref=ref(:,2:end); kp=find(sym~='z');

sym=sym(1:kp(length(kp))); % Remove predictably insignificant if kkk

sym=[sym (ref(:,kkk))']; % This is also a problem check end

seq=seqt; % Copy back the original sequence T=T/2

seq=seq.*[abs(seq)>=T==1];

seq(sig_coeff)=0; % Only those coefficients not yet found to be significant are scanned &

% coefficients previously found to be significant are made zero

num_of_passes=num_of_passes+1 drawnow

ztr=[]; % Check out end % End of while symlen=[sym 0]; save symlen %finalexamd

（3）EZW解码程序

clc format long load indices load lena

orig=x(50:50+127,50:50+127); clear x enc=symlen; T=512;

xm=mapping_128; xm=xm(:); x=xm;

dec(128*128)=0; kk=1; ii=1; c=enc(ii); root_indices=[0]; sig_ele=[0]; req=[]; aaa=0; while T>=2

% Check Dominant pass first root_indices=[0]; sig_ele=[0]; kk=1;

while (c~='0' & c~='1' & ii

if(sum([root_indices==kk])==0 & sum([find(dec~=0)==kk])==0) c=enc(ii); if c=='p'

dec(kk)=mean([T 2*T]); sig_ele=[sig_ele kk]; kk=kk+1; elseif c=='n'

dec(kk)=-1*mean([T 2*T]); sig_ele=[sig_ele kk]; kk=kk+1; elseif c=='z' Tc(kk)=0; kk=kk+1; elseif c=='r'

if kk>512 & kk<1025

temp=indices(indices(kk,1:end),1:end); root_indices=[root_indices indices(kk,1:end)

(temp(:))'];

kk=kk+1;

elseif kk>1024 & kk<4097

root_indices=[root_indices indices(kk,1:end)]; kk=kk+1; end end ii=ii+1; else kk=kk+1; end end

req=[req sig_ele(2:end)];

% Req contains all the coefficient indices in order which are found significant

if(ii0),s=1; else s=0; end

if enc(ii-1)=='1'

dec(req(rr))=abs(dec(req(rr)))+T/4; elseif enc(ii-1)=='0'

dec(req(rr))=abs(dec(req(rr)))-T/4; end

if s==0,dec(req(rr))=-dec(req(rr));end ii=ii+1; end end ii=ii-1;

disp('-------------------------------------------------------'); T

rec=round(reshape(dec(xm(:)),128,128)); xr=synth2d(rec,3,3); figure,imshow(mat2gray(xr)) PSNR(xr,orig),ii aaa=aaa+length(sig_ele)

disp('-------------------------------------------------------'); T=T/2; drawnow c=enc(ii); end

实验六 JPEG编解码

一、实验题目：

JPEG编解码器

二、实验目的：

熟悉JPEG基本系统的图像编解码方法。

三、实验内容：

编程实现近似的JPEG基本系统压缩编、解码。

四、预备知识：

（1）熟悉图像预测编码原理。 （2）熟悉图像DCT变换编码原理。 （3）熟悉熵编码原理。

（4）熟悉在MATLAB环境下对图像文件的I/O操作。

五、实验原理：

JPEG是连续色调静止图像压缩的国际标准，由于其高压缩比，使得它广泛应

用于计算机和通信等领域，例如图像压缩、多媒体通信、图像数据库等。

JPEG在压缩与解压缩的处理过程中，一般采用无损和有损两种方式。无损方式压缩比较低，有损方式则能提供很高的压缩比。 JPEG具有适中的计算复杂性，从而使得压缩算法既可以用软件实现，也可以用硬件实现。 JPEG标准中实际定义了3种编码系统：第一种是基于DCT的有损编码基本系统（顺序模式），可用于绝大多数压缩应用场合；第二种是用于高压缩比、高精度或渐进重建应用的扩展编码系统（递增模式、分层编码）；第三种是用于无失真应用场合的无损系统（DPCM）。

（1）基于 DCT的有损编码基本系统

最常用的JPEG编码是基于DCT变换的顺序型模式，又称为基本系统。其JPEG压缩编解码算法框图如图3.4.1所示。

图3.4.1 JPEG压缩编码－解码算法框图

JPEG压缩编码算法的主要计算步骤如下： ①正向离散余弦变换（ FDCT）。 ②量化(quantization)。

③使用差分脉冲编码调制（DPCM）对直流系数（DC）进行编码。 ④使用行程长度编码（RLE）对交流系数（AC）进行编码。 ⑤熵编码。

下面分别对各个计算步骤作简要介绍：

①正向离散余弦变换

在基于DCT的顺序模式中，8×8样本块以从左到右，从上到下（每次一个块行）的顺序对图像进行扫描和编码，这一过程持续到整个图像的所有数据都处理完为止。

在图3.4.1（a）的编码过程中，输入分量的样本被分成 8×8大小的数据块，并且用正向DCT（FDCT）把每个块变成64个DCT系数值，系数值中有一个是DC系数，其它63个是AC系数。其非零元素主要集中于某一个区域，通常在左上角，而右下角大部分是零。

②量化

量化是对经过FDCT变换后的频率系数进行量化。其目的是减少非“0”系数的幅度以及增加“0”值系数的数目。

由于DCT系数包含了空间频率信息，可充分利用人眼对不同频率其敏感程度不同这一特性来来选择量化表中的元素值的大小，对视觉重要的系数采用细量化（量化步长较小），如低频系数被细化；对高频系数则可采用粗量化（量化步长较大）。

JPEG给出了两个量化表，如图3.4.2所示。图（a）示出了根据视觉加权函数得到的亮度分量量化矩阵，即亮度量化表，图（b）是其色度分量量化矩阵，即色度量化表。这两个量化表是JPEG所推荐的。用户也可以设置自己的量化表，量化表必须作为JPEG编码器的输入。

量化表中的每个元素值为1～255之间的整数值，其规定了它所对应DCT系数的量化步长。

（a）亮度量化表 （b）色度量化表

图3.4.2 JPEG推荐的量化表

量化表中的某个对应值用于对相应的系数进行量化处理。量化所得数值为

零可舍去。系数量化是个十分重要的过程，它是造成DCT编解码信息损失（或失真）的根源。在JPEG中采用均匀量化器，量化定义为对64个DCT系数除以其量化步长，四舍五入取整，即：

Q（u,v）=Integer Round（F（u,v）/S（u,v））

式中，Q（u,v）为量化的系数幅度，S（u,v）为量化步长，它是量化表的元素， 量化表元素通常随DCT系数的位置和彩色分量的不同而取不同的值，量化表的尺寸为8×8与DCT 64个系数一一对应。

反量化是量化的逆过程，即

Q’=Q（u,v）×S（u,v）

③直流系数（DC）的编码

8×8 图像块经过DCT变换后，其频率系数低频分量都集中在左上角，其中F(0,0)（即第一行第一列元素）代表了直流（DC）系数，即8×8子块的平均值。直流系数有两个特点：一是系数的数值比较大，二是相邻8×8 图像块的DC系数值变化不大。根据这个特点，JPEG算法使用了差分脉冲编码调制（DPCM）技术，对相邻图像块之间的DC系数的差值单独进行编码，以提高压缩比。

④交流系数（AC）的编码

DCT变换的其他63个元素是交流（AC）系数，经量化后，其系数包含有许多“0”系数，并且许多“0”是连续的，因此使用行程长度编码（RLE）对它们进行编码。

为了保证低频分量先出现，高频分量后出现，以增加行程中连续“0”的个数，这63个元素采用了“Z”字型（Zig-Zag）的排列方法，如图3.4.3所示。

图3.4.3 量化DCT系数的编排

这63个AC系数编码的码字用两个字节来表示。如图3.4.4所示。

图3.4.4 行程编码

第一个字节的高4位表示两个非零值之间连续零的个数（行程 RunLength），低4位表示编码下一个非零值所需要的位数（尺寸Size）。第二个字节表示的是下一个非零系数的实际值。

⑤熵编码

在上面的预处理工作中，我们得到了DC码字和AC码字，为了进一步提高压缩比，需要对其再进行熵编码。熵编码过程中可使用两种编码方法： Huffman编码和自适应二进制算术编码（Adaptive Binary Arithmetic Coding）。在基本系统中推荐了两组Huffman码表，一组用于亮度信号Y，另一组用于色度信号U、V，每一组表又包括两张表，一个用于DC分量，一个用于AC分量。

编码的过程分成两步： ⅰ熵编码的中间格式表示

上述DC码字和AC码字一起作为描述符送到熵编码器进行熵编码。 对于AC系数，有两个符号。符号1为行程和尺寸，即上面的（RunLengh，Size）。（0，0）和（15，0）是两个比较特殊的情况。（0，0）表示块结束标志（EOB），

（15，0）表示ZRL，当行程长度超过15时，用增加ZRL的个数来解决，所以最多有三个ZRL（3×16＋15＝63）。符号2为幅度值（Amplitude）。

对于DC系数，也有两个符号。符号1为尺寸（Size），符号2为幅度值（Amplitude）。

ⅱ熵编码

对于AC系数，符号1和符号2分别进行编码。零行程长度超过15个时，有一个符号（15，0），块结束时只有一个符号（0，0）。

对符号1进行Huffman 编码，对符号2进行变长整数VLI编码。举例来说：Size=6时，Amplitude的范围是－63~－32，以及32~63。绝对值相同符号相反

的码字之间为反码关系，所以AC系数为32的码字为100000，33的码字为100001，－32的码字为011111，－33的码字为011110。符号2的码字紧接于符号1的码字之后。

对于DC系数，符号1和符号2也要分别进行编码。符号1用Huffman 编码，符号2用变长整数VLI编码。下面举例说明上述的编码过程。

下面为8×8的亮度（Y）图像子块经过DCT变换和量化后的系数：

15 0 -1 0 0 0 0 0 -2 -1 0 0 0 0 0 0 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

可见量化后只有左上角的几个点（低频分量）不为零，这样采用行程编码就很有效。

第一步，将AC和DC系数表示为熵编码的中间格式。先看DC系数。假设前一个8×8子块DC系数的量化值为12，则本块DC系数于它的差为3，根据下表：

Size Amplitude 0 0 1 2 3 4 5 6 7 8 9 10

-1，1 -3，-2，2，3 -7~ -4，4~7 -15~ -8，8~15 -31~ -16，16~31 -63~ -32，32~63 -127~ -64，64~127 -255~ -128，128~255 -511~ -256，256~511 -1023~ -512，512~1023

本文来源：https://www.bwwdw.com/article/67k6.html