实验九：图像编码实验
一、实验目的
1、理解图像质量评价指标及含义； 2、掌握变换图像编码原理与编程方法； 3、掌握近似JEPG图像编码原理与编程方法；
5、利用MATLAB语言自己实现常用图像编码算法，并对恢复重构的图像进行质量评价。
二、实验设备：
1、PC机或笔记本电脑，Windows操作系统； 2、安装Matlab软件，准备实验相关图像或M文件；
3、按照如下实验内容与步骤，编写相应Matlab程序，完成实验任务。 三、实验原理
JPEG编码原理与步骤（彩色图像编码，后续实验中只做灰度图像编码，其原理是相同的）： 1、 数据分块
对图像进行编码前，将每个分量图像分割成不重叠的8×8像素块， 每一个8×8像素块称为一个数据单元（DU）。在彩色图像中， JPEG分别压缩图像的每个彩色分量。虽然JPEG可以压缩通常的红绿蓝分量，但在YCbCr空间的压缩效果会更好。这是因为人眼对色彩的变化不如对亮度的变化敏感，因而对色彩的编码可以比对亮度的编码粗糙些，这主要体现在不同的采样频率和量化精度上。因此，编码前一般先将图像从RGB空间转换到YCbCr空间，再把各分量图像分割成8×8数据块。
在对图像采样时，可以采用不同的采样频率，这种技术称为二次采样。由于亮度比色彩更重要，因而对Y分量的采样频率可高于对Cb、Cr的采样频率，这样有利于节省存储空间。常用的采样方案有YUV422和YUV411。把采样频率最低的分量图像中一个DU所对应的像区上覆盖的所有各分量上的DU按顺序编组为一个最小编码单元（MCU）。对灰度图像而言，只有一个Y分量，MCU就是一个数据单元。而对彩色图像而言，以4:1:1的采样方案为例，则一个MCU由4个Y分量的DU、1个Cb分量的DU和1个Cr分量的DU组成。 2、 DCT处理
图像数据块分割后，即以MCU为单位顺序将DU进行二维离散余弦变换。对以无符号数表示的具有P位精度的输入数据， 在DCT前要减去2P-1，转换成有符号数，而在IDCT后，应加上2P-1，转换成无符号数。对每个8×8的数据块DU进行DCT后， 得到的64个系数代表了该图像块的频率成分，其中低频分量集中在左上角， 高频分量分布在右下角。系数矩阵左上角的叫做直流（DC）系数，它代表了该数据块的平均值，其余63个叫交流（AC）系数。 3、系数量化
在DCT处理中得到的64个系数中，低频分量包含了图像亮度等主要信息。在从空间域到频域的变换中，图像中的缓慢变化比快速变化更易引起人眼的注意， 所以在重建图像时，低频分量的重要性高于高频分量。因而在编码时可以忽略高频分量， 从而达到压缩的目的，这也是量化的根据和目的。
在JPEG标准中，用具有64个独立元素的量化表来规定DCT域中相应的64个系数的量化精度，使得对某个系数的具体量化阶取决于人眼对该频率分量的视觉敏感程度。理论上，对不同的空间分辨率、数据精度等情况，应该有不同的量化表。不过，一般采用下图所示的量化表，可取得较好的视觉效果。
43
之所以用两张量化表，是因为Y分量比Cb和Cr更重要些，因而对Y采用细量化，而对Cb和Cr采用粗量化。量化就是用DCT变换后的系数除以量化表中相对应的量化阶后四舍五入取整。由于量化表中左上角的值较小，而右下角的值较大， 因而起到了保持低频分量、 抑制高频分量的作用。
17182447999999991821266699999999242656999999999947669999999999999999999999999999999999999999999999999999999999999999999999999999
1611101624405161121214192658605514131624405769561417222951878062182237566810910377243555648110411392496478871031211201017292959811210010399
亮度量化表 色度量化表
4、 Z型扫描
DCT系数量化后，构成一个稀疏矩阵，用Z（Zigzag）形扫描将其变成一维数列，将有利于熵编码。Z形扫描的顺序如图：
5、 DC系数编码
DC系数反映了一个8×8数据块的平均亮度，一般与相邻块有较大的相关性。JPEG对DC系数作差分编码，即用前一数据块的同一分量的DC系数作为当前块的预测值，再对当前块的实际值与预测值的差值作哈夫曼编码。?
若DC系数的动态范围为-1024～+1024，则差值的动态范围为-2047～+2047。如果为每个差值赋予一个码字， 则码表过于庞大。因此，JPEG对码表进行了简化，采用“前缀码（SSSS）+尾码”来表示。前缀码指明了尾码的有效位数B， 可以根据DIFF从表10-8中查出前缀码对应的哈夫曼编码。尾码的取值取决于DC系数的差值和前缀码。如果DC系数的差值DIFF大于等于0，则尾码的码字为DIFF的B位原码；否则，取DIFF的B位反码。 6、 AC系数编码
经Z形排列后的AC系数，更有可能出现连续0组成的字符串， 从而对其进行行程编码将有利于压缩数据。JPEG将一个非零DC系数及其前面的0行程长度（连续0的个数）的组合称为一个事件。将每个事件编码表示为“NNNN/SSSS+尾码”， 其中， NNNN为0行程的长度，SSSS表示尾码的有效位数B（即当前非0系数所占的比特数），如果非零AC系数大于等于0， 则尾码的码字为该系数的B位原码， 否则， 取该系数的B位反码。?
44
由于只用4位表示0行程的长度，故在JPEG编码中，最大0行程只能等于15。当0行程长度大于16时，需要将其分开多次编码， 即对前面的每16个0以“F/0”表示，对剩余的继续编码。
四、练习内容
启动电脑，打开MATLAB编程编程环境，尝试建立以下M文件，且认真阅读相关代码，分析其功能，按“第五部分的要求”，完成相关实验：
1、行程编码(RLE）
clear all;
I=imread('girl.bmp'); %读入图像数据 [zipped,info]=RLEncode(I); %调用RLE进行编码 unzipped=RLEdecode(zipped,info) %调用解码程序进行解码
%显示原始图像和经编码解码后的图像,显示压缩比,并计算均方根误差得erms=0,表示%RLE是无失真编码
subplot(121);imshow(I);
subplot(122);imshow(unzipped); erms=compare(I(:),unzipped(:)) cr=info.ratio
whos I unzipped zipped
%说明：编码子程序
function [zipped,info]=RLEncode(vector)
[m,n]=size(vector); vector=vector(:)';
vector=uint8(vector(:)); L=length(vector);
c=vector(1);e(1,1)=c;e(1,2)=0; %e(:,1)存放灰度,e(:,2)存放行程 t1=1 for j=1:L
if((vector(j)==c))
e(t1,2)=double(e(t1,2))+1; else
c=vector(j); t1=t1+1; e(t1,1)=c; e(t1,2)=1; end end
zipped=e; info.rows=m; info.cols=n; [m,n]=size(e);
info.ratio=m*n/(info.rows*info.cols);
%说明：解码子程序
function unzipped=RLEdecode(zip,info)
zip=uint8(zip); [m,n]=size(zip); unzipped=[]; for i=1:m
section=repmat(zip(i,1),1,double(zip(i,2))); unzipped=[unzipped section]; end
unzipped=reshape(unzipped,info.rows,info.cols); %%程序结束
下面程序是对cameraman.tif二值化后的图像用行程编码：
45
I=imread('cameraman.tif');
I=im2bw(I,0.4); %图像二值化 [zipped,info]=RLEncode(I);
unzipped=RLEdecode(zipped,info); I2=logical(unzipped) erms=compare(I,I2) cr=info.ratio
2、变换编码
例1考虑一幅大小512×512，灰度级256的标准图像Lena，用FFT实现图像数据的压缩。首先将图像分割成(512/8)2个8×8子图像，对每个子图像进行FFT，这样每个子图像有64个傅里叶变换系数。按照每个系数的方差来排序，由于图像是实值的，其64个复系数只有一半有差别的。舍去小的变换系数，就可以实现数据压缩。这里我们保留32个系数，实现2:1的数据压缩，然后进行逆变换的程序。其Matlab程序如下：
% 设置压缩比cr
cr=0.5; % cr=0.5为2:1；cr=0.125为8:1 % 读入并显示原始图像
I1=imread('lena.bmp'); % 图像的大小为512×512
I1=double(I1)/255; % 图像为256级灰度图像，对图像进行归一化操作 figure(1); % 显示原始图像 imshow(I1);
% 对图像进行FFT
fftcoe=blkproc(I1,[8 8],'fft2(x)'); %将图像分割成8×8的子图像进行FFT coevar=im2col(fftcoe,[8 8],'distinct'); %将变换系数矩阵重新排列 coe=coevar;
[y,ind]=sort(coevar); [m,n]=size(coevar);
snum=64-64*cr; % 根据压缩比确定要变0的系数个数 % 舍去不重要的系数 for i=1:n
coe(ind(1:snum),i)=0; % 将最小的Snum个变换系数清0 end
B2=col2im(coe,[8 8],[512 512],'distinct'); % 重新排列系数矩阵
%对子图像块进行FFT逆变换获得各个子图像的恢复图像，并显示压缩图像
I2=blkproc(B2,[8 8],'ifft2(x)'); % 对截取后的变换系数进行FFT逆变换 figure(2); % 显示压缩后图像 imshow(I2);
%计算均方根误差erms
e=double(I1)-double(I2); [m,n]=size(e);
erms=sqrt(sum(e(:).^2)/(m*n))
例2 考虑例1相同的图像。首先将图像分割成(512/8)2个8×8子图像，对每个子图像进行DCT，这样每个子图像有64个变换系数，舍去50%小的变换系数，即保留32个系数，进行2:1的压缩。其Matlab程序如下：
% 设置压缩比cr
cr=0.5; % cr=0.5为2:1；cr=0.125为8:1 % 读入并显示原始图像
initialimage=imread('lena.bmp');
initialimage=double(initialimage)/255; figure(1); imshow(initialimage); % 对图像进行DCT t=dctmtx(8);
dctcoe=blkproc(initialimage,[8 8],'P1*x*P2',t,t'); coevar=im2col(dctcoe,[8 8],'distinct'); coe=coevar;
46
[y,ind]=sort(coevar); [m,n]=size(coevar); % 舍去不重要的系数 snum=64-64*cr; for i=1:n
coe(ind(1:snum),i)=0; end
b2=col2im(coe,[8 8],[512 512],'distinct'); % 对截取后的变换系数进行DCT逆变换 i2=blkproc(b2,[8 8],'P1*x*P2',t',t); % 显示压缩图像 figure(2); imshow(i2)
%计算均方根误差
e=double(initialimage)-double(i2); [m,n]=size(e);
erms=sqrt(sum(e(:).^2)/(m*n))
例3仍使用例1中的图像。首先将图像分割成(512/8)2个8×8子图像，对每个子图像进行哈达玛变换，这样每个子图像有64个变换系数，舍去50%小的变换系数，即保留32个系数，进行2:1的压缩。其Matlab程序如下：
% 设置压缩比cr
cr=0.5; % cr=0.5为2:1；cr=0.125为8:1 % 读入并显示原始图像
I1=imread('lena.bmp'); % 图像的大小为512×512
I1=double(I1)/255; % 图像为256级灰度图像，对图像进行归一化操作 figure(1); % 显示原始图像 imshow(I1);
% 对图像进行哈达玛变换
T=hadamard(8); %用于产生8×8的哈达玛矩阵 htCoe=blkproc(I1,[8 8],'P1*x*P2',T,T); CoeVar=im2col(htCoe,[8 8],'distinct'); Coe=CoeVar;
[Y,Ind]=sort(CoeVar); [m,n]=size(CoeVar); % 舍去不重要的系数 Snum=64-64*cr; for i=1:n
Coe(Ind(1:Snum),i)=0; end
B2=col2im(Coe,[8 8],[512 512],'distinct'); % 对截取后的变换系数进行哈达玛逆变换 I2=blkproc(B2,[8 8],'P1*x*P2',T, T); % 显示压缩图像 I2=I2./(8*8); figure(2); imshow(I2);
%计算均方根误差
e=double(I1)-double(I2); [m,n]=size(e);
erms=sqrt(sum(e(:).^2)/(m*n))
附1：算术编码、解码 function y=arithmetic()
format long e; clear all;
format long e;
47