利用Matlab进行图像的编码与压缩(仿照jpeg)

本文章包含以下内容：

1．选择matlab自带图像cameraman.tif和westconcordorthophoto.png为编码压缩的测试图像。

2．根据设计要求选择图像编码和压缩方法；

3、设计编码与压缩的算法；

4、显示原始图像和压缩重构图像；

5、计算压缩的性能指标：压缩率和保真度；

6、画出压缩率和保真度的关系曲线。

本文章中有：

        图像读取，余弦变换，亮度量化，Z字形重排，游程编码，参数文件保存，

        参数文件读取，逆游程编码，逆Z字形重排，逆量化压缩，生成压缩比，保真度。

仿照JPEG压缩代码如下：

（可以放到同一个文件里，R2009b以上可以用，不然的话要改函数）

如果觉着数据不好，可以更改Q_New函数中的值。

因为用到了blockproc函数，在R2009b版本以上可以直接用，否则需要改函数。

因为为了画曲线，进行了100+次图像压缩，运行时间较长，想要时间短点，可以更改36行的for循环。

用到了行程编码，所以如果图像复杂，可能会越压越大。

PadPartialBlocks作用是裁剪时，如果不足8*8，则填充成8*8的矩阵

function U()
    clear;
    clc;
    img = imread('cameraman.tif');  % 读取文件
%     img = imread('westconcordorthophoto.png');  % 读取文件
    Q = 30;  % 品质（1~100），100最好
    Q_ = Q_New(Q);
    % blkproc对图像进行分块处理,分成8*8的子块，并用dct2()处理
    DCT_c = blkproc(img,[8,8],'dct2');
    figure();     % 开一个窗口
    subplot(2,3,1);imshow(img);title('原图');   %显示原图
    % 显示余弦变换结果
    subplot(2,3,2);imshow(log(abs(DCT_c)+1),[0,10]);title('余弦变换结果');
    T = Y_Table();  % 用亮度量化表,用于量化压缩
    DCT_c = blkproc(DCT_c,[8,8],'round(x./(P1*P2))',T,Q_);   % 量化压缩
    % 显示量化压缩结果
    subplot(2,3,3);imshow(log(abs(DCT_c)+1),[0,10]);title('量化压缩结果');
    Z = blockproc(DCT_c,[8,8],@Z_open,'PadPartialBlocks',true); % Z字形重排（展开）Z_open
    S = size(Z);    % 记录展开后大小，方便以后复原
    Z = Z(:);       % 展开成一维
    Z = forecast(Z);    % 预测编码
    Z = RLC(Z);     % 游程编码
%     save('图像.mat','Q','S','Z');   % 将大小，游程编码后结果，保存到文件中
%     load('图像.mat');           % 读取保存的变量
    Z = iRLC(Z);    % 游程解码
    Z = iforecast(Z);    % 预测编码
    Z = reshape(Z,S(1),S(2));      % 重组
    Z = blockproc(Z,[64,1],@Z_merge);           % 逆Z字重排（合并）
    DCT_c = blkproc(Z,[8,8],'x.*P1*P2',T,Q_);    % 逆量化压缩
    % 显示逆量化压缩结果
    subplot(2,3,6);imshow(log(abs(DCT_c)+1),[0,10]);title('逆量化压缩结果');
    % 显示逆变换图像
    subplot(2,3,4);imshow(uint8(blkproc(DCT_c,[8,8],'idct2')));
    title('逆变换图像');
    s = size(img);  % 得到原图像的大小，用于计算压缩率
    s_ = s(1)*s(2);
    x_=[];
    y_=[];
    for i = 1:1:100
        [S,Z] = IC(img,i);  % 压缩
        s_Z = size(Z);      % 压缩后大小
        x_(length(x_(:))+1) = s_/(s_Z(1)*s_Z(2) + 3);   % 压缩比，+3是要储存游程编码大小和Q
        img_=iIC(i,S,Z);
        img_=img_(1:s(1),1:s(2));
        y_(length(y_(:))+1) = psnr(uint8(img_),img);   % 信噪比，客观保真度
        % figure,imshow(uint8(img_));
    end
    % 绘制压缩比-保真度
    subplot(2,3,5);plot(x_,y_);xlabel('压缩比');ylabel('客观保真度');
    title('压缩比-保真度');
end
 
% 由1~100,映射到实际用的Q上
function Q_ = Q_New(Q)
% 品质（1~100），100最好
    if(Q>100)   % 大于一百的和小于1的都拉回来
        Q = 100;
    end
    if(Q<1)
        Q = 1;
    end
    Q_ = 3.03-0.03*Q;    % 真正用的Q是这个，你可以改
end
 
% 图像压缩编码函数，输入图像，品质（1~100），100最好得到压缩后的数据
function  [S,Z] = IC(img,Q)
    Q_ = Q_New(Q);
    % blkproc对图像进行分块处理,分成8*8的子块，并用dct2()处理
    DCT_c = blkproc(img,[8,8],'dct2');
    T = Y_Table();  % 用亮度量化表压缩
    DCT_c = blkproc(DCT_c,[8,8],'round(x./(P1*P2))',T,Q_);
    Z = blockproc(DCT_c,[8,8],@Z_open,'PadPartialBlocks',true); % Z字形重排（展开）Z_open
    S = size(Z);    % 记录展开后大小，方便以后复原
    Z = Z(:);       % 展开成一维
%     Z = forecast(Z);    % 预测编码
    Z = RLC(Z);     % 游程编码
end
 
% 图像解压函数,输入游程编码结果,Z字重排前大小,输出解码后图像
function img = iIC(Q,S,Z)
    Q_ = Q_New(Q);
    T = Y_Table();  % 亮度量化表
    Z = iRLC(Z);    % 游程解码
%     Z = iforecast(Z);    % 预测解码
    Z = reshape(Z,S(1),S(2));  % 重组
    Z = blockproc(Z,[64,1],@Z_merge);   % 逆Z字重排（合并）
    Z = blkproc(Z,[8,8],'x.*P1*P2',T,Q_);
    img = blkproc(Z,[8,8],'idct2');
end
 
% 掏出一个亮度量化表
function T = Y_Table()
T =[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];
end
 
% 掏出Z字形重排的地址(8*8)
function T = Z_address()
T =[ 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 ];
end
 
% Z字展开（矩阵变一维）
function Z_out = Z_open(Z_input)
    Z_out = zeros(64,1);
    Z_out(Z_address()) = Z_input.data;
end
 
% Z字合并（一维变矩阵）
function Z_out = Z_merge(Z_input)
    Z_out = Z_input.data(Z_address());
end
 
% 游程编码
function Y = RLC(X)
    % 手动处理第一个
    Y(1) = X(1);
    Y(2) = 1;
    c = 1;          % 计数器，记录当前信号在Y的位置
    C = 2;          % 当前信号的游程，关于我为什么不用c+1，这样好看。
    s = max(size(X));    % 看看长度，一会有用 
    for i = (2:s)   % 对于每个数据
        if(X(i)~=Y(c))  % 如果和当前的不一样
            c = C + 1;  % 信号位置改变
            C = c + 1;  % 游程位置改变
            Y(c) = X(i);    % 记录信号
            Y(C) = 1;       % 计数1
        else
            Y(C) = Y(C) + 1;    % 游程加一
        end
    end
end
 
% 游程解码
function Y = iRLC(X)
    c = 1;          % 计数器，记录Y该到哪里了
    s = max(size(X));    % 看看长度，一会有用 
    for i = (1:2:s) % 记录格式为 信息 游程
        for j = (1:X(i+1))  % 每个信息加游程个数个
            Y(c) = X(i);    % 信息录入
            c = c + 1;      % 计数+1
        end
    end
end

% 预测函数
function Y = forecast(X)
    if(length(X)<3)
        Y = X;
        return 
    end
    s1 = X(1);
    s2 = X(2);
    Y(1)=X(1);
    Y(2)=X(2);
    for i = (3:length(X))   % 线性预测
        Y(i) = X(i)-(s1+s2)/2;
        s1 = s2;
        s2 = X(i);
    end
%     Y(50:90)
end

% 预测解码
function X = iforecast(Y)
    if(length(Y)<3)
        X = Y;
        return 
    end
    s1 = Y(1);
    s2 = Y(2);
    X(1)=Y(1);
    X(2)=Y(2);
    for i = (3:length(Y))   % 线性预测
        X(i) = Y(i)+(s1+s2)/2;
        s1 = s2;
        s2 = X(i);
    end
end

结果示例：