JPEG图像压缩性能分析(附Matlab完整代码)
注:
1、本实验不生成、不存储具体编码,只计算编码长度、PSNR和压缩比等,即计算JPEG性能。
2、本文只提供代码。如果需要完整的实现过程,压缩包下载地址:
一、实验准备
实验环境:Matlab R2017b
实验准备:lenaXING.mat含:
①codelength.mat:霍夫曼编码码长矩阵(AC、DC亮度编码表),个人整理;
②lena512.mat from https://www.ece.rice.edu/~wakin/images/;
③lighttable.mat:JPEG标准亮度量化表;
④zigzag.mat:快速Z型排序行向量。
二、实验流程
为了简化步骤,我们省略了第一步,而直接对源图像进行处理。
1、图像分割
以8x8为最小单元分割,可分割成4096个方块,从上往下存。得到32768x8的矩阵。
2、DCT变换
对这4096个方阵分别进行DCT变换,得到4096个变换方阵,从上往下存。仍是32768x8的矩阵。
3、量化
对这4096个方阵分别根据JPEG亮度标准量化表进行量化,从上往下存。仍是32768x8的矩阵。
由于四舍五入,许多高频分量已经被舍弃,减小了视觉冗余。
4、系数编码和熵编码
进一步地,我们要对图像进行系数编码和熵编码。
先用一个例子说明以下过程:a. 对每一个方阵采用ZIG-ZAG排序,以增加图像中的连0个数。得到4096x64的矩阵。
b. 对第一列的DC分量进行DPCM编码。第一个记为0,后面的都只记录与前者的差值。
c. 然后,我们根据VLI编码表,对这些数值进行分组。
d. 我们将组号记录下来。例如:(2,3)意思是:DPCM编码为3,且位于第二组(如表所示)。
DC组号(实际上就是编码长度)记录在4096x1的矩阵里。
e. 根据组号和DC值,计算出VLI和霍夫曼编码码长和。
f. 我们对其余63列AC分量进行RLC编码。
首先,我们要得到这些AC分量的中间格式。假设-30前面有1个0,并且属于VLI编码表中的第五组,则其中间格式为: (1,5),-30。
注意,若连0超过15,则记作(15,0)。
g. 根据中间格式,计算总编码长度。
(这个表是我根据AC亮度霍夫曼表整理的,可能有错误。如果需要使用该表,请务必校正)
三、实验结果
压缩倍数约为16.8。
峰值信噪比为35.8。
四、实验代码
function LenaJPEGbyXING % 数字图像处理:lena512图像的JPEG压缩 % by XING, Beijing clear;close all;clc; %% Image Segmentation % lena512: 512*512 % Block: 32768*8 load('lenaXING.mat','lena512'); Block=[]; for numi=1:64 %逐行取方阵 m=(numi-1)*8+1; %每块行的开头 for numj=1:64 %逐列取方阵 n=(numj-1)*8+1; %每块列的开头 Block=[Block; lena512(m:m+7,n:n+7)]; end end %% DCT % Block: 32768*8 % FBlock: 32768*8 for num=1:4096 start=(num-1)*8+1; FBlock(start:start+7,:)=dct2(Block(start:start+7,:)); end %% Quantization % FBlock: 32768*8 % QBlock: 32768*8 load('lenaXING.mat','lighttable'); for num=1:4096 start=(num-1)*8+1; QBlock(start:start+7,:)=round(FBlock(start:start+7,:)./lighttable); end %% Inverse Quantization % QBlock: 32768*8 % reFBlock: 32768*8 for num=1:4096 start=(num-1)*8+1; reFBlock(start:start+7,:)=QBlock(start:start+7,:).*lighttable; end %% IDCT for num=1:4096 start=(num-1)*8+1; Block(start:start+7,:)=idct2(reFBlock(start:start+7,:)); end %% Image Reconstrucion relena512=[]; for numi=1:64 m=(numi-1)*512+1; % 分成64个512*8阵列,每个阵列有64个8*8方阵 A=[]; for numj=1:64 n=(numj-1)*8; A=[A Block(m+n:m+n+7,:)]; end relena512=[relena512; A]; end %% JPEG & JPEG2000 Figure figure(); subplot(1,2,1); imshow(lena512./256); xlabel('Origin'); subplot(1,2,2); imshow(relena512./256); xlabel('JPEG'); suptitle('Origin vs. JPEG by XING'); figure(); subplot(1,2,1); imshow(relena512./256); xlabel('JPEG'); subplot(1,2,2); lena2k = imread('lena512.bmp'); imwrite(lena2k,'lena_16.8.j2k','CompressionRatio',16.8); imshow('lena_16.8.j2k') xlabel('JPEG2000'); suptitle('JPEG vs. JPEG2000 by XING'); %% PSNR delta=lena512-relena512; delta=delta.^2; MSE=sum(delta(:))/512/512; PSNR=10*log10(255^2/MSE); disp(['PSNR_JPEG: ',num2str(PSNR)]); lena16_8=imread('lena_16.8.j2k'); delta=lena2k-lena16_8; delta=delta.^2; MSE=sum(delta(:))/512/512; PSNR=10*log10(255^2/MSE); disp(['PSNR_JPEG2000: ',num2str(PSNR)]); %% CODE % 以下只计算编码长度,不实际存储编码。 %% ZIG-ZAG % QBlock: 32768*8 % QLine: 4096*64 QLine=[]; load('lenaXING.mat','zigzag'); zigzag = zigzag + 1; % 下标加1,从0开始 for num=1:4096 start=(num-1)*8+1; A=reshape(QBlock(start:start+7,:),1,64);% 变成行向量 QLine=[QLine;A(zigzag)]; end %% DPCM for DC % QLine: 4096*64 % VLIDC: 4096*1 % 对第一列进行DPCM编码,第一个值记为DC,并赋0 DC=QLine(1,1);%保留备用 sumcode=0;%计算编码长度 QLine(1,1)=0; for num=4096:-1:2 QLine(num,1)=QLine(num,1)-QLine(num-1,1); end VLIDC=ones(4096,1);% VLI分组 for num=1:4096 temp=abs(QLine(num,1));%用绝对值判断组别 if temp==0 VLIDC(num)=0; else for k=1:7%经测试,第一列最大值为80,前7组够用 if (temp>=2^(k-1)) && (temp<2^k) VLIDC(num)=k; break; end end end end for num=1:4096 %先根据DC亮度huffman表计算sumcode if (VLIDC(num)<=5) && (VLIDC(num)>=0) sumcode=sumcode+3; elseif VLIDC(num)==6 sumcode=sumcode+4; else sumcode=sumcode+5; end %再根据VLI表计算sumcode sumcode=sumcode+VLIDC(num); end %DC计算结果为24096,比8*4096=32768小得多。 %% RLC for AC % QLine: 4096*64 % 经测试,后63列最大值为58,VLI前6组够用。 eob=max(QLine(:))+1; %设一个超大值作为每一行结束符 for numn=1:4096 %放eob for numm=64:-1:2 if QLine(numn,numm)~=0 QLine(numn,numm+1)=eob; break; end if (numm==2)%没找到 QLine(numn,2)=eob; end end end test=QLine'; [col,~]=find(test==eob);%我们只要eob列位置 validAC=col-1; %每一行保留的AC数据量,含EOB maxcz=0; for numn=1:4096 %逐行计算并加至sumcode cz=[];%记录前0数 VLIAC=[];%记录组号 count=0;%记录连0数 for numm=2:1+validAC(numn) if QLine(numn,numm)==eob cz=[cz 0]; VLIAC=[VLIAC 0]; elseif QLine(numn,numm)==0 count=count+1; else %遇到非0值 if count>15 %遇到连0大于15的 cz=[cz 15]; count=0; VLIAC=[VLIAC 0]; continue; end cz=[cz count]; count=0; temp=abs(QLine(numn,numm));%用绝对值判断组别 for k=1:6%经测试,后63列最大值为58,前6组够用 if (temp>=2^(k-1)) && (temp<2^k) VLIAC=[VLIAC k]; break; end end end end%该行cz和VLIAC已定,开始计算 sumcode=sumcode+4;%EOB对应1010,就是4bit czlen=length(cz)-1; %czlen不包括EOB load('lenaXING.mat','codelength'); for k=1:czlen if VLIAC(k)==0 sumcode=sumcode+11; else sumcode=sumcode+codelength(cz(k)+1,VLIAC(k)); end end end %sumcode计算结果为124555。 %% Compression Rate OB=512*512*8; CR=OB/sumcode; PD=sumcode/512/512; disp(['Original Bit: ',num2str(OB),' bit']); disp(['Compressed Bit: ',num2str(sumcode),' bit']); disp(['Compression Ratios: ',num2str(CR)]); disp(['Pixel Depth: ',num2str(PD),' bpp']); disp(' ——Calculated by XING'); end