Transformation For Image Compression二维信号

这是二维信号处理(2-D)课程的Term,简单做了几种二维信号变换算法的对比 。
图像压缩算法就是尽可能去除信号里的冗余信息 。从可视化的角度来说冗余信息就是不影响视觉效果的信息,这种信息一般来说是图像的边缘(edge)信息 。想象一下绘画过程 , 图像可以分为色块()和边缘(edge),如果去掉边缘信息只保留色块 , 仍然能看出原图的样子(就相当于画画不勾边).如果去掉色块信息,图像会变得模糊(不能等同于画画只勾边不上色,去掉色块信息以后把一些空间相关信息也去掉了 , 不仅仅只是去掉颜色) , 参照下图 。
从二维信号的角度,信息属于低频信息,edge信息属于高频信息 。上图就是使用高斯低通滤波器对图像进行处理(PS里的高斯模糊就是这个算法),图像变得模糊,出现“失真” 。图像压缩应尽可能保留低频信息,过滤高频信息 。
对二维信号进行变换,目的是使得图像的能量变得集中,能量稀疏的部分经过舍弃就变为0,0越多编码方便 , 存储编码比存储原图的矩阵节省空间 。
注意,图像压缩不仅仅只是对二维信号进行变换处理 , 还需要经过量化和编码,这里不讨论这两个步骤 , 主要讨论变换处理 。
算法介绍
类似于一维信号的傅里叶变换,二维信号(一般以图像形式呈现)也有相应的二维傅里叶变换2-D(2D DFT). 和一维信号类似,2D DFT就是把信号分解成许多正弦函数和余弦函数之和 。这些余弦函数和正弦函数被称为基函数(Basis ) 。如果把基函数换成其他的函数 , 那么也可以有对应形式的分解变换 。下面介绍五种分解变换的方式,分别是2-D,2-DSine
,2-DWalsh–,2-D和。
2-D(DCT)
DCT 变换的基函数就是函数,变换以后对信号进行对称延展,所以频谱具有连续性 。这也是著名的图像压缩算法JPEG(就是那个常用图片格式 , 其实这是一种压缩算法的简写)所使用的分解变换算法 。DCT变换后的信号,能量全部集中在频谱的左上低频区,方便编码和保留 。
JPEG算法中,图像不是以全图为单位进行处理,而是先把图像分成了*的小块,对每个小块分别进行变换 。我们也模仿这种形式,先分块,后变换 , 下图为变换后的效果,可以看到 , 能量都集中在每个小块的左上角(白色明亮区域代表高能量区域) 。
2-DSine(DST)
类似DCT,DST就是把基函数换成了Sin 。DST有很多不同的定义,我选了其中一种,这种定义下的DST变换频谱不连续,低频区不集中,能量的分布也被打散了,如图 。
2-DWalsh–(DWHT)
DWHT的基函数不是三角函数,而是方波函数(但不一定是正方形 , 方波的宽度可以随意),方波的值无非+1和-1,所以这个方法计算起来很方便 。结果类似于DCT,如下图.
2-D(DWT)
之前学一维变换的时候,小波变换我就没理解 。学二维的时候居然能懂一点了,小波也是波 , 但是区别于无穷的、周期性的正弦波余弦波,它是长度有限的,波形不一定很规则 。二维小波变换就是把图像分解成许许多多长度有限的小波,小波的种类还有很多种,这里只用了其中一种 , 叫 。
二维小波变换的操作相当于两次滤波和降采样,第一次对行(row)进行高通或低通滤波 , 然后降采样,行的长度变为原来的一半,第二次对列进行高通或低通滤波 , 然后降采样,列的长度又变为原来的一半 。所得的子图会有四种情况:经历两次高通滤波、经历先高通后低通滤波、经历先低通后高通滤波、经历两次低通滤波 。
下图是对全图进行小波变换,没有分块 。DWT已经被应用在新的图像压缩算法中,即算法. 该算法不像JPEG,需要先将图像分成许许多多小块,它可以直接对原图进行变换,或者对比较大的分块(比如1/4图大?。┙斜浠淮?。
变换后所得到的数字(其实是不同小波分量的系数),可以分为四组:
一组构成(四张的左上),这组经历了两次高通滤波,去除了行和列上的高频信息,可以出,如前面所说 , 去除高频信息对图像的影响最小 。
一组构成 (右上), 这组对行(row)进行低通滤波,对列进行高通滤波,所以保留的是垂直方向上的重要信息,显现出来的就是垂直方向的低频高能量分量 。
一组构成 (左下), 这组对行(row)进行高通滤波 , 对列进行低通滤波,和的情况正好相反,保留水平方向上的重要低频信息 。
最后一组是  , 这组把水平和垂直方向上的低频信息都过滤了,保留最多的是对角线上的信息 。
下图是分块处理的结果 , 可以更清楚地看到被保留的水平、垂直、对角线低频高能量信息 。能量的分布就是沿着水平、垂直、对角线方向 。
2-D(PCT)
PCT这个名称不是最常见的,比较常见的名称是-Loève ,或者。
这个变换比较特殊,不是基于将原信号投射到基函数上,是基于找出矩阵的特征值,“正交化” 。
一维中 , 这个变换就是把一个non-zero mean、向量变成一个zero-mean、的向量
在二维图像处理中,这个变换可以去除图像不同的“层”之间的相关性(这就是不同层之间的冗余信息 , 层的一个常见形式是RGB, 彩色图像由R\G\B三个色层组合构成) 。
下图是分解后的结果,把一个RGB图分解成三个低相关的正交矩阵(一维信号是向量,二维信号就是矩阵) 。
从每个层(即图上的band,即通道)的图像中,可以得到的一个信息是:右上P1通道的信息含量最多,包括了space、 and edge,P2通道的信息量少于P1,但是仍然有一些关键信息,P3通道之所以看起来全黑,是因为它几乎不含什么信息 。
彩图压缩算法的实现(附代码及解释)
**读入图像,将数据形式转变为以便于后续计算 。**是调用的函数包(后面会展示函数包里的函数):
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就是变换后的二维频谱图 。
()是自带的函数,可以实现对矩阵的分块处理 。
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变换的主程序:
().()是调用了自定义函数,后面会详解 。这里对得到的系数进行了四舍五入处理(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);endfor 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三个输出就是变换后得到的三个通道,E是的矩阵
[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}); %irreversibleend
编码以及存储编码数据:
编码也会极大影响压缩率,这里就用了最常见的赫夫曼编码和零游程编码(zero run- )的结合,然后转换成可存储的01形式 。()函数是从网上找的一个压缩存储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);endv_com=v_com(:)';[code1,dict1] = imgprocess().encoding(v_com);[ zipped, pad ] = save_binary( code1 );fid=fopen("compressed_file.txt","w");fwrite(fid,zipped,"uint8");
解码并转变成矩阵
从这步开始进行解压缩 。注意省去了读取压缩文件的步骤,直接,然后整理成矩阵 。DWT 和PCT的代码和其它三种变换不同 。
ticSig{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,峰值信噪比,用来衡量重构图的质量 。
figureimshow(uint8(I))figureimshow(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彩图并分成三个通道分别处理
某些函数不能对三维矩阵进行处理,所以必须分成三个矩阵一个个处理 。
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 indexM = log2(N)+1;binHadIdx = fliplr(dec2bin(HadIdx,M))-'0'; % Bit reversing of the binary indexbinSeqIdx = zeros(N,M-1);% Pre-allocate memoryfor k = M:-1:2% Binary sequency index binSeqIdx(:,k) = xor(binHadIdx(:,k),binHadIdx(:,k-1));endSeqIdx = binSeqIdx*pow2((M-1:-1:0)');% Binary to integer sequency indexwalshMatrix = hadamardMatrix(SeqIdx+1,:) ;% 1-based indexingc = 1./(walshMatrix)/sqrt(N);end
图像分块代码
写了一个代码,把图像分成正方形或者长方形小块,不满足指定大小的部分(一般就是边缘区域)则保留成其原本的大?。ū热纾?9的图像分成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;endendblocks=blocks-1; %最后一次执行后仍然加一,-1后才是实际的block数量end
子图拼接复原
上面那个函数的反过程,把分出来的小块拼合成原来的图

Transformation For Image Compression二维信号

文章插图
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];endendend
分块进行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, 可以组成处理后的图,可用于encodeA_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 vectorif isnumeric(x)x = x(:);x = x(~isnan(x));xid = [];else[x,xid] = grp2idx(x);x = x(~isnan(x));endx = sort(x(:));m = length(x);x1 = diff(x);x1(end + 1) = 1;x1 = find(x1);CumFreq = x1/m;value = http://www.kingceram.com/post/x(x1);x1 = [0; x1];Freq1 = diff(x1);Freq2 = Freq1/m;end
编码
结合了零游程编码的赫夫曼编码,用了自带的赫夫曼编码函数
function [code,dict] = encoding(Sig) x=Sig;%zero run-length codingy=[];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;endend[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 runc=[];i=1;while i<=length(r)if r(i)==0 && i
压缩存储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型
结果示例
自己用三张不同的图试了一下 , 结果作为参考 。每一张图中包括了原图和压缩并复原的图 。
参考文献
装模作样列上参考文献,毕竟为了做这个还是读了一点东西的 。
[1] M. , Y. , M. L. , A. , and T. A. , “ evalu-ation of DWTto DCT forimage,”ofand, vol. 6, no. 4, pp. 9–15, Apr. 2014, doi: 10.5815/.2014.04.02.
A. H. M. J. I. , T. A. , and K. , “Anfor Color Image -sion of JPEG and PNGUsing DCT and DWT,” 2014on -and, pp. 129–133, Nov. 2014, doi: 10.1109/cicn.2014.40.
[2] S. A. Joshi, H. N. , A. , V. Patil, and A. P. , “of 2D-DCT based JPEG,” 2023onand Sus-(), pp. 1–6, Jun. 2023, doi: 10.1109/.2023..
[3] M.and S. , Newfor JPEG . 2020. doi: 10.1109/.2020..
[4] S.and D. , Eds., “Newfor JPEG ,” 33rdon Com-puter , May 2017, doi: 10.48550/arXiv.1705.03531.
[5] Furht, E. Akar, and W. A. , “ JPEG image ,” inin, 2018, pp. 35–50. doi: 10.1007/978-3-319-96634-2_6.
[6] Y. L. , Y. , Y. , and R. , “ oflossy image,” 2021on,Com- and Smart(ICSES), Sep. 2021, doi: 10.1109/.2021..
[7] U.and C. H. Chen,Image : Anwith . 2009. []. :
[8] V. Tyagi,image . 2018. doi: 10.1201/.
[9] E.and R. M. ,. 1983. []. -ble:
[10] W.and M. J. Burge,ofimage :.&Media, 2013.
[11] R. C.and R. E. Woods,Image . 2002.
[12] M. Sonka, V. , and R. Boyle, Image ,and. , 2014.
[13] M.and G. ,and imageusing . John Wiley & Sons, 2010.
[14] S. E. ,imageand :withand . CRC Press, 2017.
[15] J. W. Woods,, image, and videoand .Press, 2011.
【Transformation For Image Compression二维信号】![在这里插入图片描述]( # =)