图像DCT 变换编码与压缩
一、实验目的：
（1）掌握离散余弦变换DCT 的实现方法，了解DCT 的幅度分布特性，从而加深对DCT 变换的认识。
（2）掌握图像DCT 变换编码的实现方法，从而加深对变换编码压缩图像原理的理解。
二、实验内容：
编程实现图像DCT 变换编码。
三、实验原理：
变换编码是在变换域进行图像压缩的一种技术。
图1显示了一个典型的变换
编码系统。
压缩图像
输入图像N×N
图1 变换编码系统
在变换编码系统中，如果正变换采用DCT 变换就称为DCT 变换编码系统。
DCT 用于把一幅图像映射为一组变换系数，然后对系数进行量化和编码。
对于大多数的正常图像来说，多数系数具有较小的数值且可以被粗略地量化（或者完全抛弃），而产生的图像失真较小。
DCT 是仅次于K-L 变换的次最佳正交变换，且以获得广泛应用，并成为许多图像编码国际标准的核心。
离散余弦变换的变换核为余弦函数，计算速度快，有利于图像压缩和其他处理。
对于N ×N 的数字图像，二维DCT 变换的正反变换可表示为：
11
0011
00
(21)(21)(,)()()(,)cos
cos 222
(21)(21)(,)()()(,)cos
cos 22N N x y N N u v x u y v F u v c u c v f x y N N x u y v f x y c u c v F u v N
N N
ππ
ππ
--==--==++=++=
∑∑∑∑（1）
其中，
1/00()()1,1,2,...,1
u v c u c v u v N ⎧==⎪==⎨
=-⎪⎩或 MATLAB 图像处理工具箱实现离散余弦变换有两种方法：
（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 的矩阵进行计算的函数，P i是传递给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显示了采用JPEG标准化矩阵进行DCT变换编码的结果。
图2 DCT变换编码
四、实验步骤：
DCT变换编码流程如下:
步骤1：设置JPEG标准化数组;
步骤2：求8×8快的DCT变换矩阵;
步骤3: 计算8×8快的DCT变换；
步骤4：对DCT系数量化和反量化；
步骤5：求反量化系数的逆DCT变换；
步骤6：重新显示重建图像、误差图像和误差图像的直方图。
量化时可采用JPEG标准推荐的归一化数组，如表1所示。
表1 JPEG标准化数组
程序如下：
x=imread('lena','bmp');
[M,N]=size(x);%得到原始举矩阵大小
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 10;
72 92 95 98 112 100 103 99];
I=x;
x=double(x);%变成双精度
t=dctmtx(8);%得到DCT变换矩阵
y1=blkproc(x,[8 8],'P1*x*P2',t,t');%进行DCT变换得到变换矩阵
y2=blkproc(y1,[8 8],'round(x./P1)',m);%量化
y3=blkproc(y2,[8 8],'x.*P1',m);%反量化
y=blkproc(y3,[8 8],'P1*x*P2',t',t);%反DCT变化IDCT
subplot(2,2,1);imshow(I);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(:),256);%生成直方图数据
subplot(2,2,4);bar(k,h,'k');title('误差图像直方图');
五、思考题目：
（1）观察图像8×8子块的DCT系数的分布，并分析其特点。
答：经DCT变换以后，系数大多数集中在左上角（即低频分量），其余系数大多很小或为零。
（2）将量化步长分别增大为初始值的2倍、4倍、8倍后再进行DCT变换编码，显示不同量化步长条件下的重建图像、误差图像以及误差图像的直方图。
分析重建图像质量和量化步长的关系。
答：由以上结果对比可看出随着量化补偿的加大，图像误差变大，失真越来越严重。