JPEG图像压缩 matlab实现

注：

1、该博客已经停止维护了，相关文件也早已丢失了，各位非常抱歉！

2、本实验不生成、不存储具体编码，只计算编码长度、PSNR 和压缩比等，即计算 JPEG 性能。

3、本文只提供代码。如果需要完整的实现过程，压缩包下载地址：

https://download.csdn.net/download/weixin_41730407/10371917

4、该程序还使用 imwrite 实现了相同放大倍数下 JPEG2000 的压缩，并与 JPEG 压缩指标进行了对比。

实验环境：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 LenaJPEG
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