Transformation For Image Compression二维信号变换算法在图像压缩上的应用比较

这是二维信号处理（2-D Signal Processing）课程的Term Project，简单做了几种二维信号变换算法的对比。

图像压缩算法就是尽可能去除信号里的冗余信息。从可视化的角度来说冗余信息就是不影响视觉效果的信息，这种信息一般来说是图像的边缘（edge）信息。想象一下绘画过程，图像可以分为色块（region）和边缘（edge），如果去掉边缘信息只保留色块，仍然能看出原图的样子（就相当于画画不勾边）.如果去掉色块信息，图像会变得模糊（不能等同于画画只勾边不上色，去掉色块信息以后把一些空间相关信息也去掉了，不仅仅只是去掉颜色），参照下图。

从二维信号的角度，region信息属于低频信息，edge信息属于高频信息。上图就是使用高斯低通滤波器对图像进行处理（PS里的高斯模糊就是这个算法），图像变得模糊，出现“失真”。图像压缩应尽可能保留低频信息，过滤高频信息。
对二维信号进行变换，目的是使得图像的能量变得集中，能量稀疏的部分经过舍弃就变为0，0越多编码方便，存储编码比存储原图的矩阵节省空间。
注意，图像压缩不仅仅只是对二维信号进行变换处理，还需要经过量化和编码，这里不讨论这两个步骤，主要讨论变换处理。

算法介绍Algorithm Introduction

类似于一维信号的傅里叶变换，二维信号（一般以图像形式呈现）也有相应的二维傅里叶变换2-D Discrete Fourier Transform（2D DFT）. 和一维信号类似，2D DFT就是把信号分解成许多正弦函数和余弦函数之和。这些余弦函数和正弦函数被称为基函数（Basis Function）。如果把基函数换成其他的函数，那么也可以有对应形式的分解变换。下面介绍五种分解变换的方式，分别是2-D Discrete Cosine Transform，2-D Discrete Sine
Transform，2-D Discrete Walsh–Hadamard Transform，2-D Discrete Wavelet Transform 和 Principal Components Transform。

2-D Discrete Cosine Transform （DCT）

DCT 变换的基函数就是Cosine函数，变换以后对信号进行对称延展，所以频谱具有连续性。这也是著名的图像压缩算法JPEG（就是那个常用图片格式，其实这是一种压缩算法的简写）所使用的分解变换算法。DCT变换后的信号，能量全部集中在频谱的左上低频区，方便编码和保留。
JPEG算法中，图像不是以全图为单位进行处理，而是先把图像分成了8pixel*8pixel的小块，对每个小块分别进行变换。我们也模仿这种形式，先分块，后变换，下图为变换后的效果，可以看到，能量都集中在每个小块的左上角（白色明亮区域代表高能量区域）。

2-D Discrete Sine Transform （DST）

类似DCT，DST就是把基函数换成了Sin。DST有很多不同的定义，我选了其中一种，这种定义下的DST变换频谱不连续，低频区不集中，能量的分布也被打散了，如图。

2-D Discrete Walsh–Hadamard Transform （DWHT）

DWHT的基函数不是三角函数，而是方波函数（但不一定是正方形，方波的宽度可以随意），方波的值无非+1和-1，所以这个方法计算起来很方便。结果类似于DCT，如下图.

2-D Discrete Wavelet Transform （DWT）

之前学一维变换的时候，小波变换我就没理解。学二维的时候居然能懂一点了，小波也是波，但是区别于无穷的、周期性的正弦波余弦波，它是长度有限的，波形不一定很规则。二维小波变换就是把图像分解成许许多多长度有限的小波，小波的种类还有很多种，这里只用了其中一种，叫Daubechies。
二维小波变换的操作相当于两次滤波和降采样，第一次对行（row）进行高通或低通滤波，然后降采样，行的长度变为原来的一半，第二次对列进行高通或低通滤波，然后降采样，列的长度又变为原来的一半。所得的子图会有四种情况：经历两次高通滤波、经历先高通后低通滤波、经历先低通后高通滤波、经历两次低通滤波。
下图是对全图进行小波变换，没有分块。DWT已经被应用在新的图像压缩算法中，即JPEG2000算法. 该算法不像JPEG，需要先将图像分成许许多多小块，它可以直接对原图进行变换，或者对比较大的分块（比如1/4图大小）进行变换处理。

变换后所得到的数字（其实是不同小波分量的系数），可以分为四组：
一组构成approximation（四张subimage的左上），这组经历了两次高通滤波，去除了行和列上的高频信息，可以出，如前面所说，去除高频信息对图像的影响最小。
一组构成Vertical Detail(右上), 这组对行（row）进行低通滤波，对列进行高通滤波，所以保留的是垂直方向上的重要信息，显现出来的就是垂直方向的低频高能量分量。
一组构成Horizontal Detail(左下), 这组对行（row）进行高通滤波，对列进行低通滤波，和Vertical的情况正好相反，保留水平方向上的重要低频信息。
最后一组是Diagonal Detail，这组把水平和垂直方向上的低频信息都过滤了，保留最多的是对角线上的信息。
下图是分块处理的结果，可以更清楚地看到被保留的水平、垂直、对角线低频高能量信息。能量的分布就是沿着水平、垂直、对角线方向。

2-D Principal Components Transform （PCT）

PCT这个名称不是最常见的，比较常见的名称是Karhunen-Loève Transform，Principal Components Transform或者Eigenvalue Transform。
这个变换比较特殊，不是基于将原信号投射到基函数上，是基于找出矩阵的特征值，“正交化”。
一维中，这个变换就是把一个non-zero mean、correlated向量变成一个zero-mean、uncorrelated的向量
在二维图像处理中，这个变换可以去除图像不同的“层”之间的相关性（这就是不同层之间的冗余信息，层的一个常见形式是RGB, 彩色图像由R\G\B三个色层组合构成）。
下图是分解后的结果，把一个RGB图分解成三个低相关的正交矩阵（一维信号是向量，二维信号就是矩阵）。

从每个层（即图上的band，即通道）的图像中，可以得到的一个信息是：右上P1通道的信息含量最多，包括了space、region and edge，P2通道的信息量少于P1，但是仍然有一些关键信息，P3通道之所以看起来全黑，是因为它几乎不含什么信息。

彩图压缩算法的Matlab实现（附代码及解释）

**读入图像，将数据形式转变为double以便于后续计算。**imgprocess是调用的函数包（后面会展示函数包里的函数）：

clear;
imgprocess; img=imread("building1.tif");
[imR, imG, imB, imRGB] = imgprocess().imread3(img);
%I = im2double(I(:,:,1));
I = double(imRGB);

实现DCT\DST\DWHT的主程序：
这三个算法的实现形式很接近，可以共用大部分代码，所以写在一起。变换的计算方法是矩阵乘法，具体自己学习。三个变换只是生成了不同的矩阵，后面的乘法计算相同。im_T就是变换后的二维频谱图。
blockproc（）是matlab自带的函数，可以实现对矩阵的分块处理。

N=16;
for i=1:3im_T=I(:,:,i);T = dctmtx(N);%T = imgprocess().dwhtmtx(N) ;%T = imgprocess().dstmtx(N) ;formula = @(block_struct) T * block_struct.data * T'; im_T = blockproc(im_T,[N N],formula);B(:,:,i) = im_T;
%imshow(B)
end

分块计算DWHT变换的主程序：
imgprocess().blockdwt（）是调用了自定义函数，后面会详解。这里对得到的系数进行了四舍五入处理（round（）），是为了方便后面的计算。这对压缩率和保真度有影响，如果想要无损压缩，就不能改变变换所得系数。
A V H D 是四个系数矩阵，R G B三个通道各有四个系数矩阵

n=432;
m=640;
for i=1:3[A_img{i}, H_img{i}, V_img{i}, D_img{i}]=imgprocess().blockdwt(I(:,:,i),n,m);
end
for i=1:3A{i} = round(A_img{i}(:))';H{i} = round(H_img{i}(:))';V{i} = round(V_img{i}(:))';D{i} = round(D_img{i}(:))';
end
[A_img H_img, V_img, D_img]=dwt2(I,'db2','mode','per');

PCT变换主程序：
同样进行了四舍五入处理。
P, C, T三个输出就是变换后得到的三个uncorrelated通道，E是eigenvector的矩阵

 [P,C,T,E,e] = imgprocess().pct(I);Sig{1}=round([P(:)]');Sig{2}=round([C(:)]');

四舍五入处理数据（仅针对DST DCT DWHT）
为了方便编码，要把矩阵处理成一维行向量，前面也有类似的操作。

v1=B(:,:,1);
v2=B(:,:,2);
v3=B(:,:,3);
v{1} = [v1(:)]';
v{2} = [v2(:)]';
v{3} = [v3(:)]';
for i=1:3Sig{i}=round(v{i}); %irreversible
end

编码以及存储编码数据：
编码也会极大影响压缩率，这里就用了最常见的赫夫曼编码和零游程编码（zero run-length coding）的结合，然后转换成可存储的01形式。save_binary（）函数是从网上找的一个压缩存储01编码的函数，不是自己写的。
一般的彩图要存储三个通道的编码信息，PCT可以只存储两个通道的编码信息，这样就减少了占用的空间，所以PCT是非常高效的压缩变换。

 [code1,dict1] = imgprocess().encoding(Sig{1});[code2,dict2] = imgprocess().encoding(Sig{2});
[code3,dict3] = imgprocess().encoding(Sig{3});
code=[code1 code2 code3];
%code=[code1 code2]; % 2 PCT coefficients matrix
[ zipped, pad ] = save_binary( code );
fid=fopen("compressed_file.txt","w");
fwrite(fid,zipped,"uint8");

DWT的整合和编码
DWT的编码程序要单独写，因为变换后得到了很多矩阵，这些矩阵要全部压成一个行向量进行编码。

A_img=round(A_img);
H_img=round(H_img);V_img=round(V_img);D_img=round(D_img);
for i=1:3%v_com=cat(1,v_com,A{i},H{i},V{i},D{i});v_com=cat(1,A_img,H_img,V_img,D_img);
end
v_com=v_com(:)';
[code1,dict1] = imgprocess().encoding(v_com);
[ zipped, pad ] = save_binary( code1 );
fid=fopen("compressed_file.txt","w");
fwrite(fid,zipped,"uint8");

解码并转变成矩阵
从这步开始进行解压缩。注意省去了读取压缩文件的步骤，直接decode，然后整理成矩阵。DWT 和PCT的代码和其它三种变换不同。

tic
Sig{1} = imgprocess().decoding(code1,dict1);
Sig{2} = imgprocess().decoding(code2,dict2);
Sig{3} = imgprocess().decoding(code3,dict3);
B=[];
for i=1:3B(:,:,i)=reshape(Sig{i},size(I(:,:,1)));
end
%____________________DWT__________________________%
% Sig{1} = imgprocess().decoding(code1,dict1);
% for i=1:3
%     A_img{i}=reshape(A{i},size(D_img{3}));
%     H_img{i}=reshape(H{i},size(D_img{3}));
%     V_img{i}=reshape(V{i},size(D_img{3}));
%     D_img{i}=reshape(D{i},size(D_img{3}));
% end
%__________PCT________________%
% Sig{1} = imgprocess().decoding(code1,dict1);
% Sig{2} = imgprocess().decoding(code2,dict2);
% P=reshape(Sig{1},size(I(:,:,1)));
% C=reshape(Sig{2},size(I(:,:,1)));

反变换重构图像
这里有个不完美的地方，就是没有研究出分块DWT反变换（前面有做分块的DWT变换），所以只能对全图进行反变换。

%________________DCT-DST-DWHT_________________________________%
invt = @(block_struct) T' * block_struct.data * T;
I2=[];
for i=1:3I2(:,:,i) = blockproc(B(:,:,i),[N N],invt);
end
%________________DWT_________________________________%% re_img= idwt2(A_img, H_img, V_img, D_img,"db2"); 
%________________PCT_________________________________%
% re_img=imgprocess().ipct(P, C, E);

展示重构的图
这里还计算了PSNR，峰值信噪比，用来衡量重构图的质量。

figure
imshow(uint8(I))
figure
imshow(uint8(I2));title("Decompressed image - with 4 \times 4 blocks")
%imshow(uint8(re_img));title("Decompressed image- with no tiled")
PSNR=psnr(uint8(I2),uint8(I))
%I = imresize(I,size(re_img(:,:,1)));
%PSNR=psnr(uint8(re_img),uint8(I))
toc

下面是一些用到的自定义函数。

读取RGB彩图并分成三个通道分别处理
Matlab某些函数不能对三维矩阵进行处理，所以必须分成三个矩阵一个个处理。

function [imR, imG, imB, imRGB]=imread3(img)imR = img(:,:,1);imG = img(:,:,2);imB = img(:,:,3);imRGB = img;
end

构造DST变换中用到的矩阵

function c = dstmtx(n)
validateattributes(n,{'double'},{'integer' 'scalar'},mfilename,'n',1);[cc,rr] = meshgrid(0:n-1);c = sqrt(2 / n) * sin(pi * (cc + 1/2) .* (rr+1) / ( n));
c(1,:) = c(1,:) / sqrt(2);end

构造DWHT变换中用到的矩阵

function c = dwhtmtx(N)
hadamardMatrix = hadamard(N);
HadIdx = 0:N-1;                          % Hadamard index
M = log2(N)+1;
binHadIdx = fliplr(dec2bin(HadIdx,M))-'0'; % Bit reversing of the binary index
binSeqIdx = zeros(N,M-1);                  % Pre-allocate memory
for k = M:-1:2% Binary sequency index binSeqIdx(:,k) = xor(binHadIdx(:,k),binHadIdx(:,k-1));
end
SeqIdx = binSeqIdx*pow2((M-1:-1:0)');    % Binary to integer sequency index
walshMatrix = hadamardMatrix(SeqIdx+1,:) ;% 1-based indexing
c = 1./(walshMatrix)/sqrt(N);
end

图像分块代码
写了一个代码，把图像分成正方形或者长方形小块，不满足指定大小的部分（一般就是边缘区域）则保留成其原本的大小（比如：99的图像分成44，那就得到 4个44 和 4个14 还有一个1*1）。

function  [block,blocks,row_blk_num,col_blk_num]=decompose(img,n,m)
%n*m is the size of block
[row, col] = size(img);
row_blk_num = ceil(row / n);   
col_blk_num = ceil(col / m);   block = cell(1,row_blk_num*col_blk_num);
blocks = 1;
for i = 1:row_blk_numfor j = 1:col_blk_num%disp(blocks);%只为了显示计数，可以删掉block{blocks} = img((i - 1) * n + 1 : ((i*n > row)*end + (i*n <= row)*i*n), ...(j - 1) * m + 1 : ((j*m > col)*end + (j*m <= col)*j*m), :);blocks = blocks + 1;end
endblocks=blocks-1; %最后一次执行后仍然加一，-1后才是实际的block数量
end

子图拼接复原
上面那个函数的反过程，把分出来的小块拼合成原来的图

function re_img=compose(row_blk_num,col_blk_num,block)%处理后的子图按顺序和大小复原%dwt2()
%wcodemat(:,255,'mat',1);validateattributes(block,{'cell'},{'nonempty'},mfilename,'n',1);blocks = 1;
re_img = [];
temp = [];
for i = 1:row_blk_numtemp = [];for j = 1:col_blk_numtemp2 = block{blocks};temp = [temp,temp2];blocks = blocks + 1;endif(i==1)re_img = temp;endif(i~=1)re_img=[re_img;temp];  end
end
end

分块进行DWT

function [A_img, H_img, V_img, D_img, row_blk_num,col_blk_num]=blockdwt(img,n,m)[block,blocks,row_blk_num,col_blk_num]=imgprocess().decompose(img,n,m);
Vcell=cell(1,blocks); %预分配空间
Hcell=cell(1,blocks); 
Dcell=cell(1,blocks); 
Acell=cell(1,blocks); for i=1:blocks[A, H, V, D] =  dwt2(block{i},'db1','mode','per'); %transformAmtx = wcodemat(A,255,'mat',1);Vmtx = wcodemat(V,255,'mat',1); %convert to matirx 1~225Hmtx = wcodemat(H,255,'mat',1);Dmtx = wcodemat(D,255,'mat',1);Acell{i} = Amtx;Vcell{i} = Vmtx;Hcell{i} = Hmtx;Dcell{i} = Dmtx;
end%coefficient, 可以组成处理后的图，可用于encode
A_img=imgprocess().compose(row_blk_num,col_blk_num,Acell);
V_img=imgprocess().compose(row_blk_num,col_blk_num,Vcell);
H_img=imgprocess().compose(row_blk_num,col_blk_num,Hcell);
D_img=imgprocess().compose(row_blk_num,col_blk_num,Dcell);end

PCT正变换

function [P,C,T,E,e] = pct(varargin)
%----------  Input/Output Argument Check ----------------------------if nargin<1,error('Too few arguments for pct_cvip');elseif nargin>1,error('Too many arguments for pct_cvip');end;    if nargout<4,error('Too few output arguments for pct_cvip');elseif nargout>5,error('Too many output arguments for pct_cvip');end;   
%--------- RGB Image Input Check -----------------------------------if size(varargin{1},3)~=3error('Invalid Image Input: Requires Color(RGB) Image');end;    
%--------- Data Type Check and Conversion ---------------------------    if ~isa(varargin{1},'double')varargin{1}=double(varargin{1});end;    
%---------- RGB to PCT Conversion ------------------------------r=varargin{1}(:,:,1); %r bandg=varargin{1}(:,:,2); %g bandb=varargin{1}(:,:,3); %b bandp = size(r,1)*size(r,2);
%----------red mean, green mean and blue mean--------------------mr = sum(sum(r))/p;mg = sum(sum(g))/p;mb = sum(sum(b))/p;        
%Autocovariance terms of covariance matrixCrr = sum(sum((r-mr).^2))/p;Cgg = sum(sum((g-mg).^2))/p;Cbb = sum(sum((b-mb).^2))/p;
%Cross-covariance terms of covariance matrixCgr = sum(sum((g-mg).*(r-mr)))/p;Cbr = sum(sum((b-mb).*(r-mr)))/p;Cbg = sum(sum((b-mb).*(g-mg)))/p;
%Covariance matrixCOVrgb = [Crr Cgr Cbr;Cgr Cgg Cbg;Cbr Cbg Cbb];
%EigenVector and EigenValues (E and e)[E,e] = eigs(COVrgb);E = E';        
% %Components of PCTP = r*E(1,1)+g*E(1,2)+b*E(1,3);C = r*E(2,1)+g*E(2,2)+b*E(2,3);T = r*E(3,1)+g*E(3,2)+b*E(3,3);end

PCT反变换
用两个正交通道复原RGB

function re_img=ipct(P, C, E)E=inv(E);R = P*E(1,1)+C*E(1,2);G = P*E(2,1)+C*E(2,2);B = P*E(3,1)+C*E(3,2);re_img=cat(3,R,G,B);
end

统计矩阵中的元素及其出现频率
这个函数不是原创的。用于生成赫夫曼编码的字典。

function [value,Freq2] = HistRate(x)
%   HistRate(x)，displays a frequency table of the data in a vector
if isnumeric(x)x = x(:);x = x(~isnan(x));xid = [];
else[x,xid] = grp2idx(x);x = x(~isnan(x));
end
x = sort(x(:));    
m = length(x);
x1 = diff(x);    
x1(end + 1) = 1;
x1 = find(x1);
CumFreq = x1/m;
value = x(x1);
x1 = [0; x1];
Freq1 = diff(x1);
Freq2 = Freq1/m;
end

编码
结合了零游程编码的赫夫曼编码，用了Matlab自带的赫夫曼编码函数

function [code,dict] = encoding(Sig) 
x=Sig;%zero run-length coding
y=[];
c=1;%计数器
i=1;
while i<=length(x) if(x(i)==0)j=i+1;if (j<=length(x)) && (x(j)~=0) y=[y,x(i),1];i=i+1;elsewhile (j+1<=length(x))&& (x(j+1) == 0) j=j+1;endif j>length(x)y=[y,x(i),1];breakendnum_0=j-i+1;y=[y,x(i),num_0];i=i+num_0;endelsey=[y,x(i)];c=1;i=i+1;end
end
[symbols, p]= imgprocess().HistRate(y);dict = huffmandict(symbols',p');code = huffmanenco(y,dict);
end

解码

function c = decoding(code,dict)r = huffmandeco(code,dict);% Undo terminating zero run
c=[];
i=1;
while i<=length(r)if r(i)==0 && i<length(r)c=[c zeros(1,r(i+1))];i=i+2;elsec=[c r(i)];i=i+1;end
end
end

压缩存储01编码的函数
从网上找的

% 用Huffman编码字典对给定数据(int16)编码，并进一步编码将结果整理为(uint8)类型
function [ zipped, pad ] = save_binary( string )string = uint8(string);            % uint8型，1字节%补零: 便于每次取8个数的处理
len = length(string);
pad = 8-mod(len,8);     % 须存储！！
if pad>0   % 如果string长度不能被8整除，则在其右端补零直到能被8整除为止string = [string uint8(zeros(1, pad))];
end% 每8位存储为一个uint8型整数（0~255）
cols = length(string) / 8;          % 列数
string = reshape(string, 8, cols);  % 将string改编成8*cols的矩阵
weights = 2.^(0:7);                 % 1*8的矩阵
zipped = uint8(weights*double(string)); % 矩阵相乘，得1*cols的矩阵，unit8型

结果示例

自己用三张不同的图试了一下，结果作为参考。每一张图中包括了原图和压缩并复原的图。

参考文献

装模作样列上参考文献，毕竟为了做这个Project还是读了一点东西的。
[1] M. Beladgham, Y. Habchi, M. L. Abdelmouneim, A. Bassou, and T. A. Abdelmalik, “Performance evalu-ation of DWT compared to DCT for compression Biomedical image,” International Journal of Modern Education and Computer Science, vol. 6, no. 4, pp. 9–15, Apr. 2014, doi: 10.5815/ijmecs.2014.04.02.
A. H. M. J. I. Barbhuiya, T. A. Laskar, and K. Hemachandran, “An Approach for Color Image Compres-sion of JPEG and PNG Images Using DCT and DWT,” 2014 International Conference on Computa-tional Intelligence and Communication Networks, pp. 129–133, Nov. 2014, doi: 10.1109/cicn.2014.40.
[2] S. A. Joshi, H. N. Prajwal, A. Gumaste, V. Patil, and A. P. Guddimath, “Performance Analysis of 2D-DCT based JPEG Compression Algorithm,” 2023 International Conference on Applied Intelligence and Sus-tainable Computing (ICAISC), pp. 1–6, Jun. 2023, doi: 10.1109/icaisc58445.2023.10199973.
[3] M. Morsli and S. Yaichi, New Transforms for JPEG Format. 2020. doi: 10.1109/icmit47780.2020.9047014.
[4] S. Svoboda and D. Barina, Eds., “New Transforms for JPEG Format,” 33rd Spring Conference on Com-puter Graphics, May 2017, doi: 10.48550/arXiv.1705.03531.
[5] Furht, E. Akar, and W. A. Andrews, “Lossless JPEG image compression,” in SpringerBriefs in computer science, 2018, pp. 35–50. doi: 10.1007/978-3-319-96634-2_6.
[6] Y. L. Prasanna, Y. Tarakaram, Y. Mounika, and R. Subramani, “Comparison of different lossy image compression techniques,” 2021 International Conference on Innovative Computing, Intelligent Com-munication and Smart Electrical Systems (ICSES), Sep. 2021, doi: 10.1109/icses52305.2021.9633800.
[7] U. Qidwai and C. H. Chen, Digital Image Processing: An Algorithmic Approach with MATLAB. 2009. [Online]. Available: http://ci.nii.ac.jp/ncid/BB02040823
[8] V. Tyagi, Understanding digital image processing. 2018. doi: 10.1201/9781315123905.
[9] E. Dudgeon and R. M. Mersereau, Multidimensional digital signal processing. 1983. [Online]. Availa-ble: http://ci.nii.ac.jp/ncid/BA0035985X
[10] W. Burger and M. J. Burge, Principles of digital image processing: Advanced Methods. Springer Science & Business Media, 2013.
[11] R. C. Gonzalez and R. E. Woods, Digital Image Processing. 2002.
[12] M. Sonka, V. Hlavac, and R. Boyle, Image processing, analysis and machine vision. Springer, 2014.
[13] M. Charbit and G. Blanchet, Digital signal and image processing using MATLAB. John Wiley & Sons, 2010.
[14] S. E. Umbaugh, Digital image processing and analysis: Applications with MATLAB and CVIPtools. CRC Press, 2017.
[15] J. W. Woods, Multidimensional signal, image, and video processing and coding. Academic Press, 2011.