块变换编码--DCT编码

概念

块变换编码，是指将图像分成若干个块，对每一块经过数学转换后映射至另一值域后再进行编码处理，是属于变换编码的一种。因为经常与量化一起使用，属于有损数据压缩。

流程图

以下是块变换编码的流程图（有损也可以从压缩的步骤中看出来）：

Created with Raphaël 2.1.0block transform encoderInput image(N*N)Input image(N*N)Construct n*n subimagesConstruct n*n subimagesForward transformForward transformQuantizerQuantizerSymbol endoderSymbol endoderCompressed imageCompressed image

Created with Raphaël 2.1.0block transform decoderCompressed imageCompressed imageSymbol endoderSymbol endoderInverse transformInverse transformMerge n*n subimagesMerge n*n subimagesDecompressed imageDecompressed image

变换形式及原理

其中包括以下几种形式：

傅里叶变换

DCT变换

Walsh-Hadamard Transform

KL Transform

简单来说，就是通过把图像从空间域变为另外一种域，用频域作为例子，因为图像一般都具有低频的特性，所以其实转换成频域之后会有很多高频的成分是可以忽略并且不影响我们认为的“图像质量”（不过眼神好的人可能会发现图像边缘变模糊了），因此先作变换，再处理。

K–L变换的压缩效率很高，但算法实现困难；FFT变换算法实现简单，但压缩效率不是很理想。而DCT变换虽然不是最佳的的变换，但因为计算简单， 是常用的变换方法（至于为什么计算简单估计要从数学推导上入手了，这个我还没有推导，希望有大神能够指点一下）。

实现

下面用matlab（R2015a）实现DCT变换中的分别用区域编码和阈值编码方法实现图像压缩，用 8×8DCT 变换，保留 50％ 的系数(区域编码保留前 50%个系数，阈值编码保留 50%的大系数）。

clear;
%% 显示原始图像
A=imread('D:\\Picture\\1.jpg');%这里换成相应的路径即可
[Imgm,Imgn]=size(A);
Target = imcrop(A, [0,0,floor(Imgn/8)*8,floor(Imgm/8)*8]);%为了用blockproc函数
I=double(rgb2gray(Target));%将彩色图像变为灰度图像
figure(1);
imshow(uint8(I));
title('原始图像');
%% dct变换，并且显示图像
T=zeros(8,8); %变换矩阵A,也可以通过函数dctmtx(n)求得
[m,n]=size(T);
for i=0:m-1
for j=0:m-1
if i==0
a=sqrt(1/m);
else
a=sqrt(2/m);
end
T(i+1,j+1)=a*cos(pi*(j+0.5)*i/m);
end
end
%----------------------------------------------------------
% B = blockproc(A,[m n],fun)
% 输入灰度图像A，返回图像B，按照尺寸m× n滑动邻域，利用运算函数fun处理后得到结果
% [m n] ： 图像以m*n为分块单位，对图像进行处理（如8像素*8像素）
% fun： 应用此函数对分别对每个m*n分块的像素进行处理
%----------------------------------------------------------
fun = @(block_struct) T * block_struct.data * T';
I1=blockproc(I,[8 8],fun);
figure(2);
imshow(uint8(I1));
title('DCT变换');
%% 区域编码，并且显示图像
a=[];
a0=ones(1,8);
a1=zeros(1,8);
%区域编码保留前50%个系数
for i=1:8
if i<=4
a(i,:)=a0;
else
a(i,:)=a1;
end
end;
fun = @(block_struct) a.*block_struct.data;
areaCode=blockproc(I1,[8 8],fun);
figure(3);
imshow(uint8(I1));
title('对DCT区域编码');
%% 区域编码解码，并且显示图像
fun = @(block_struct) T' * block_struct.data * T;
areaCodeDecode=blockproc(areaCode,[8 8], fun);
figure(4);
imshow(uint8(areaCodeDecode));
title('对DCT区域编码后反DCT变换');
%% 阈值编码，并且显示图像
[Im,In]=size(I1);
%----------------------------------------------------------
% B=reshape(A,m,n)
% 返回m-by-n矩阵B，它的元素从A的第一列开始取到A的最后一列（一列取完再取第二列）
% 如果A没有m*n元素，得到一个错误结果
%----------------------------------------------------------
list=reshape(I1,1,Im*In);%将I1变成一个1*(Im*In)的列向量

%----------------------------------------------------------
% B=median(M)
% 每一列返回一个值,为M该列的从大到小排列的中间值.
%----------------------------------------------------------
mid=median(abs(list));%即取出g中绝对值的中间值，因为DCT变换正负并不影响大小

%阈值编码保留大的50%个系数
for i=1:m
for j=1:n
if(abs(I1(i,j))<mid)
a(i,j)=0;
end;
end;
end;
fun = @(block_struct) a.* block_struct.data;
thresholdCode=blockproc(I1,[8 8],fun);
figure(5);
imshow(uint8(thresholdCode));
title('对DCT阈值编码');

%% 阈值编码解码，并且显示图像
fun = @(block_struct) T' * block_struct.data * T;
thresholdCodeDecode=blockproc(thresholdCode,[8 8],fun);
figure(6);
imshow(uint8(thresholdCodeDecode));
title('对DCT阈值编码后反DCT变换');

到这里整个DCT变换就已经做完了，不过在这个过程中我还发现了一点有意思的东西，就是

B = blockproc(A,[m n],fun)

这个函数，当[m n]模板访问A完了之后还有剩余的情况下（比如说我用8*8的模板去访问1366*768的图片时，我们会发现无论是横向还是纵向都不能除尽，也就是说会有一些边缘的部分会不能填满8*8的模板），fun函数如果设置成矩阵相乘的形式，就像我上面所写的

fun = @(block_struct) T’ * block_struct.data * T

的话，那么matlab就会报错说维度不相同，不能处理。因为在一开始读入图片的时候用了

Target = imcrop(A, [0,0,floor(Imgn/8)*8,floor(Imgm/8)*8])

来规整图片的像素，丢弃边缘的信息，当然这里选取的8只能因为我的滑动窗是8*8的。

另外，用

blkproc(A,[m n],fun, parameter1, parameter2, …)

则没有这个问题，可是这个函数对于大型图片处理并不如blockproc方便（没有存储功能），大家可以自行斟酌。