数字图像处理 实验七：JPEG压缩编码

DIP实验7：JPEG压缩编码

实验目的

1）掌握n×n块的DCT图像变换及频谱特点。
2）熟悉JPEG基本系统的图像编解码方法。

实验原理

基于DCT有损压缩编码的基本系统的主要计算步骤
1）通过前向离散余弦变换（FDCT）减少图像数据相关性；
2）利用人眼的视觉特性对DCT系数进行量化；

附：可供参考的Matlab函数有dctmtx/dct2、blkproc、hist、bar、num2str和strcat

实验内容

1）编程实现n×n块DCT变换的图像频谱显示，块DCT系数按照Zig-Zag扫描并取部分进行图像重建，计算图像的均方根误差RMSE，显示误差图像和误差直方图。
2）对8×8块的DCT系数，采用JPEG默认的量化矩阵进行量化和反量化，计算原图像与重建图像之间的均方根误差RMSE、并显示误差图像。

输出图像排列格式如下：

原图像	8×8 DCT频谱显示	取?个DCT系数时的重建图像
差值直方图, RMSE=	Scale为?时的重建图像	量化误差图像, RMSE=

需要编写的函数：
1）分块DCT函数：
根据给定的n值，计算原图像的n×n块的DCT系数，并转换为可视频谱图以便观察

function dctfre = DctBlock(oldbuf,n,t);
% 输入
%     oldbuf         		 原始图像数据
%     n              		 分块的大小
%     t             	     DCT变换矩阵
% 输出
%     dctfre                 DCT系数矩阵
%     log(abs(dctfre)*5+1)   修正DCT系数以便频谱观察

2）ZigZag及IDCT函数：
按Zig-Zag扫描顺序，只保留64个DCT系数中的前x个系数，并对修改后的DCT系数用逆DCT变换重建图像。将x值从1开始增大，观察图像由模糊到清晰的变化情况。计算重建图像的均方根误差RMSE，显示误差图像和误差直方图。

function e = ZigIDCT(oldbuf,dctfre,x,n,t);
% 输入
%     oldbuf          原始图像数据
%     dctfre          DCT系数矩阵
%     x                每个分块中需要保留的DCT系数个数
%     n                分块的大小
%     t                DCT变换矩阵
% 输出
%     e                原图像与保留部分DCT系数后的重建图像之间的误差矩阵

3）量化和反量化函数：
根据给定的默认JPEG量化表，对每个n×n块的DCT系数进行量化和反量化，显示量化误差图像及其直方图。

function QuanIQuan(oldbuf,dctfre,Scale,n,t);
% 输入
%     oldbuf          原始图像数据
%     dctfre          DCT系数矩阵
%     Scale           乘积因子
%     n               分块的大小
%     t               DCT变换矩阵

zigzag = [1 2 6 7 15 16 28 29          
          3 5 8 14 17 27 30 43
          4 9 13 18 26 31 42 44
          10 12 19 25 32 41 45 54
          11 20 24 33 40 46 53 55
          21 23 34 39 47 52 56 61
          22 35 38 48 51 57 60 62
          36 37 49 50 58 59 63 64];
          
% Luminance quantization table
QuanTable =      
	[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]; 

---

参考代码

OI = imread('lena.bmp');
n = 8;    % DCT分块大小
x = 2;    % 每个分块中需要保留的DCT系数个数

dctfre = DctBlock(OI, n);    % 分块DCT函数
[Lerror, Lrestored] = ZigIDCT(OI, dctfre, x, n);    % ZigZag及IDCT函数
 
Scale = 2;    % 乘积因子
[Qerror, Qrestored] = QuanIQuan(OI, dctfre, Scale, n);    % 量化和反量化函数

subplot(2, 3, 1);     
imshow(OI);     
title('原图像');

dct_amend = log(abs(dctfre) * 5 + 1);    % 修正DCT系数以便频谱观察
dct_amend = Calibration(dct_amend);

subplot(2, 3, 2);       
imshow(dct_amend);      
title([num2str(n), '×', num2str(n), ' DCT频谱显示']);

subplot(2, 3, 3);
imshow(Lrestored);       
title(['取', num2str(x), '个DCT系数时的重建图像']);

error_lessen_rms = sqrt( sum(sum(Lerror.^2)) / numel(Lerror) );
% numel: 返回数组或者向量中所含元素的总数

subplot(2, 3, 4);
histogram(Lerror);
title(['差值直方图, RMSE=', num2str(error_lessen_rms)]);

subplot(2, 3, 5);
imshow(Qrestored);        
title(['Scale为', num2str(Scale), '时的重建图像']);
 
Qerror_rms = sqrt( sum(sum(Qerror.^2)) / numel(Qerror) );
subplot(2,3,6);
imshow(uint8(Calibration(Qerror)));     
title(['量化误差图像,RMSE=', num2str(Qerror_rms)]);


%% 分块DCT函数
function dctfre = DctBlock(oldbuf, n)
	% 输入
	%     oldbuf         原始图像数据
	%     n              分块的大小
	% 输出
	%     dctfre         DCT系数矩阵
	oldbuf = double(oldbuf);
	oldbuf = oldbuf - 128;    % DCT变换范围：-128~127
	D = dctmtx(n);    % 返回 n×n 的DCT变换矩阵
	% 矩阵A的DCT变换可用 D×A×D' 来计算
	% blkporc: 功能为对图像分块进行DCT变换
	% B = blkproc(A, [M N], Fun) ，A为输入图像，M*N为块大小，Fun为处理函数 
	dctfre=blkproc(oldbuf, [n, n], 'P1*x*P2', D, D'); 
	% D为变换矩阵( D, D'为要传递给函数的参数)，P1和P2为参数，代表D*x*D’ 
end


%% ZigZag及IDCT函数
function [e, Restored] = ZigIDCT(oldbuf, dctfre, x, n)
	% 输入
	%     oldbuf        原始图像数据
	%     dctfre         DCT系数矩阵
	%     x                 每个分块中需要保留的DCT系数个数
	%     n                 分块的大小
	% 输出
	%     e                 原图像与保留部分DCT系数后的重建图像之间的误差矩阵
	%     Restored    重建图像
	 
	zigzag = [1  2   6   7   15  16  28  29          
	          3  5   8   14  17  27  30  43
	          4  9   13  18  26  31  42  44
	          10 12  19  25  32  41  45  54
	          11 20  24  33  40  46  53  55
	          21 23  34  39  47  52  56  61
	          22 35  38  48  51  57  60  62
	          36 37  49  50  58  59  63  64];
	 
	% 根据当前x值(保留细数个数)得到“Z”字扫描的逻辑值, mask为logic类型
	mask=zigzag<=x;   
	% 返回 n×n 的DCT变换矩阵
	D=dctmtx (n);      
	% Restored为重建的压缩图像矩阵
	Restored = blkproc(dctfre, [n, n], 'P1*(x.*P2)*P3', D' ,mask,D);     
	Restored = Restored+128;
	Restored = uint8(Restored); 
	% 原图像与保留部分DCT系数后的重建图像之间的误差矩阵
	e = double(oldbuf)-double(Restored);     
end

%% 量化和反量化函数
function [e, Restored] = QuanIQuan(oldbuf, dctfre, Scale, n)
	% 输入
	%     oldbuf           原始图像数据
	%     dctfre           DCT系数矩阵
	%     Scale            乘积因子
	%     n                分块的大小
	% 输出
	%     e                量化误差
	%     Restored         重建图像
	
	% Luminance quantization table
	QuanTable = [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];     
	 
	Qvalue = blkproc(dctfre, [n n],'round(x./P1)', Scale*QuanTable);     
	% 量化
	IQvalue = blkproc(Qvalue, [n n], 'x.*P1', Scale*QuanTable);     % 反量化
	% 对经过量化和反量化后的矩阵进行逆DCT变换得到重建图像矩阵
	D = dctmtx(n);
	Restored = blkproc(IQvalue, [n n], 'P1*x*P2', D', D);
	Restored = Restored + 128;
	Restored = uint8(Restored); 
	% e为量化误差矩阵
	e = double(Restored) - double(oldbuf);
end

%% 标定
function cal = Calibration(img)
    img = double(img);
    [M, N] = size(img);
    fmin = min(min(img));
    fm = img - fmin * ones(M, N);
    fmmax = max(max(fm));
    fs = 255 * fm ./ fmmax;
    cal = uint8(fs);
end

实验结果