jpegprogram.m

灰度图像的jpeg压缩处理

close all;
clear;
clc
%时间：2007年12月4日
%仿真jpeg算法的图像压缩，并进行分析
%按照jpeg编码的流程，进行运算。
%说明：jpeg压缩步骤：
% 1.从输入图像中提取8×8块，即待处理图像图像进行不重叠的分割；分割的边界不够可以补零处理。
% 2.对每一8×8块分别进行DCT变换
% 3.量化处理
% 其中量化处理，之后还要对参数进行zigzag扫描。
% 4.编码处理

%读取待处理的图像
quality=1; %压缩质量控制因子                                                最后要改的地方A1－2^8=256

qua_array = [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;%JPEG标准量化矩阵
    
% zigzag_order = [1, 9, 2, 3, 10,17,25,18,11, 4, 5,12,19,26,33,41 ...
%                 34,27,20,13, 6, 7,14,21,28,35,42,49,57,50,43,36 ...
%                 29,22,15,8, 16,23,30,37,44,51,58,59,52,45,38,31 ...
%                 24,32,39,46,53,60,61,54,47,40,48,55,62,63,56,64];
            
% Zig-Zag scanning of AC coefficients,
%Z矩阵
z=[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];



 f_pic_read=imread('len.tif');                                                 
% 
  [mf,nf]=size(f_pic_read);
  
  f_pic=double(f_pic_read)-128;   %对图像进行偏置处理32768＝2^15                  最后要改的地方C128－32768
  
  mb=mf/8; nb=nf/8;  %mb，以及nb分别是行向以及列向的分块数

 %下面对图像进行2维的分块DCT变换
  C=dctmtx(8);%余弦变换矩阵C
  
  fun =@(x) C*x*C';
  f_dct2=blkproc(f_pic,[8 8],fun);
%   f_dct2=blkproc(f_pic,[8 8],'P1*f_pic*P2', C, C');
%   %对图像分块并且进行DCT变换,可能会出错，结果上式果然出错了，注释掉，重新用在上面的两句代替
 clear f_pic;
 %对图像进行块量化
 f_qua=round(blkproc(f_dct2,[8 8],'divq',qua_array)); %也可以向上面的一样处理
 
 clear f_dct2;
%对直流系数进行编码处理
if mb*nb > 1,                               % mb*nb表示最后分成的块数  %Fq对应的是经过DCT变换，以及量化过的矩阵
   f_dc=reshape(f_qua(1:8:mf,1:8:nf)',mb*nb,1); %把矩阵(这个矩阵对应的是直流别的系数)reshape成一个列向量，列向量里面的数值全部是直流向量
   f_dc_dpcm=dpcm(f_dc,1);                      %引用函数dpcm               通过这个可以求出直流系数的差分数值，取出所有的直流分量，进行差分编码,
else
   f_dc_dpcm=f_qua(1,1);                            %如load的图像分割后只有一个DCT块，的处理
end
clear f_dc;
dc_cof=[];                                    %直流系数的初始赋值，为空
for i=1:mb*nb,
   dc_cof=[dc_cof jdcenc(f_dc_dpcm(i))];      %引用函数jdcenc               对每一个直流差值进行huffman的查表编码，并以码流的形式输出
end

clear f_dc_dpcm;

%下面对交流系数进行编码处理
ac_seq=[];
for i=1:mb
  for j=1:nb
    tmp(z)=f_qua(8*(i-1)+1:8*i,8*(j-1)+1:8*j); %
    % tmp is 1 by 64
    eobi=max(find(tmp~=0)); %end of block index
                    % eob is labelled with 999    
    ac_seq=[ac_seq tmp(2:eobi) 999];
  end
end

clear f_qua;
ac_cof=jacenc(ac_seq); %交流系数的编码
clear ac_seq;
disp(['DC coefficient after Huffman coding has ' int2str(length(dc_cof)) ...
' bits']);
disp(['AC coefficient after Huffman coding has ' int2str(length(ac_cof)) ...
' bits']);

disp(['Compression Rate   ' num2str((length(dc_cof)+length(ac_cof))/(mb*nb*64)) '   Bits / pixel '])
disp(['Compression Ratio   ' num2str(8/((length(dc_cof)+length(ac_cof))/(mb*nb*64))) ' : 1'])


%以下是解码的部分                   

dec_ac=jacdec(ac_cof);            %引用函数jacdec进行交流解码
dec_dc=jdcdec(dc_cof);            %引用函数jdcdec进行直流解码

clear ac_cof;
clear dc_cof;

%解码之后的处理
Eob=find(dec_ac==999);             %找到每一个8×8DCT块的结束标志在解码流中的坐标，在AC的编码时将结束标志定义为999，因为在8bit灰度级中，
%没有999这个值，为了避免与灰度值混淆，要是灰度级达到10bit及以上，这个数据要修改

%对解码之后的AC系数，重新排列到用于图像重建的矩阵之中
kk=1;ind1=1;n=1;               
for ii=1:mb
    for jj=1:nb
        re_ac=dec_ac(ind1:Eob(n)-1);    %Eob(n)-1是在解码流中，每一个DCT块的交流有效数据的坐标
        ind1=Eob(n)+1;               %ind1的坐标是第n+1个DCT模块的第一个数据对应的坐标，Eob(n)返回的是第n个DCT块的结束标志在解完码的数据流中的位置
        n=n+1;                       %开始处理第n＋1模块
        ri(8*(ii-1)+1:8*ii,8*(jj-1)+1:8*jj)=dezz([dec_dc(kk) re_ac zeros(1,63-length(re_ac))]);
        kk=kk+1;                     %这里kk的作用是用来对直流的系数进行定位
    end
end

clear Eob;
clear dec_ac;

re_f_qua=round(blkproc(ri,[8 8],'idivq',qua_array));
clear ri;
fun1 =@(x) C'*x*C;
re_f_dct2=blkproc(re_f_qua,[8 8],fun1);
clear re_f_qua;
restore_pic=round(re_f_dct2+128);
clear re_f_dct2;
figure(2);
imshow(mat2gray(restore_pic)),title(' Reconstructed ')

% Calculate MSE , SNR
f_pic_read=double(f_pic_read);
de_matrix=(f_pic_read-restore_pic).^2;                % Doubt about formulae
clear restore_pic;
MSE=sum(de_matrix(:));
clear de_matrix;
PSNR=10*log10(255^2*mf*nf/MSE)       % Doubt about formulae