周俊敏
JPEG压缩格式是目前图像处理领域里面用得最广泛的一种图像压缩方式，它的实现主要分成四个步骤：
1．颜色模式转换及采样；
2．DCT变换（离散余弦变换）；
3．量化；
4．编码（有算术编码和霍夫曼编码两种，这里采用霍夫曼编码），用VB语言编程实现以上四个步骤，即完成了JPEG压缩过程，这里假设给定的源图像是一幅24位真彩色的BMP图像。
一、颜色转换及采样
1．颜色转换：对BMP图像中的颜色数据进行由RGB一YCbCr的转换，Y表示亮度，CbCr分别表示蓝色度和红色度。
转换公式：
Y=0.2990R+0.5870G+0.1140B
Cb=-0.1687R-0.3313G+0.5000B
Cr=0.5000R-0.4187G-0.0813B
这样转换以后就得到三个新的基色值，对这三个基色值来讲，都可以当作一
个独立的图像平面来进行压缩编码。
2．采样：颜色转换后，保留每一点的亮度值Y而色度值Cb，Cr则是每两点保留一点，在图像的行和列方向上都可执行颜色采样，这里采用的采样比是行列方向都是2:1:1，在行方向，每两点保留一点，列方向也是每两点保留一点，这样如果假设原来的CbCr矩阵大小为M*S，则经过2:1:1抽样之后成了M/2*s/2＝1/4M*S，只有原来的1/4了，图像大大缩小，如果想减小图像的失真度，则可行方向不抽样或列方向不抽样。
程序实现：
1．读取BMP图像信息，获取图像行像素和列像素数值，在BMP图像中，图像数据是以倒序存放的。亦即实际图像第一行资料存放在BMP图像数据矩阵的最后一行，依次类推，所以取资料的时候要从BMP图像数据矩阵的最后一行开始读起，把数据存放在新建数组(或称矩阵)的第一行，一直取完。BMP图像行像素和列像素的数值分别存于文件头信息的第18和22字节。用Get语句可得到Pwidth(图像宽度)和Phight(图像高度)的数值。
2．得到Pwidth*Phight后便得知BMP图像的像素点大小，而数据矩阵（三基色，RGB矩阵）的大小是Pwidth*Phight*3，因为每一个像素点都含有RGB三个数据，我们要处理的是数据矩阵而不是像素点矩阵。所以，新建数组的大小是(Pwidth*3)*Phight。
在BMP图像数据矩阵中，三原色RGB的排列顺序是这样的：
B. G. R. B. G. R...........R
B. G. R....................R
B........
. .
B. G. R....................R
接下来，把这一个矩阵(包含BGR)拆分成三个独立的B、G、R矩阵，得到三个新的矩阵(只包含B的矩阵，只包含G的矩阵，只包含R的矩阵)，简称为B矩阵、G矩阵、R矩阵(大小为Pwidth*Phight)，用程序实现拆分，只要依次取原矩阵的第1、4、7、10、13......个资料即得到B矩阵，依次读取第2、5、8、11......个数据即得到G矩阵，依次读第3、6、9、12......个资料即得到R矩阵。
得到B、G、R矩阵后再利用颜色转换公式很容易就可得到YCbCr矩阵。
Y(n)=0.114B(n)+0.587G(n)+0.299R(n)
Cb(n)=0.5B(n)-0.3313G(n)-0.1687R(n)
Cr(n)=0.0813B(n)-0.14187G(n)+0.5R(n)
(For n=1 To PW*PH)
其中PW为Pwidth的简写，PH为Phight的简写。
得Y、Cb、Cr矩阵后，先判断一下三个矩阵的行数和列数是否都被16整除，
如果不能则以最后一行或一列的资料填充到能被16整除为止，所需填充的数目
为：行方向：16-(Pwidth mod 16)
列方向：16-(Phight mod 16)
填充以后YCbCr矩阵的大小成为M*S，其中，设M=Pwidth+16-(Pwidth mod
16);S=Phight+16-(Phight mod 16)。
3、抽样：Y矩阵不动，CbCr矩阵行列方向上相邻两点只保留一点的资料值
(只要在其行方向取第1、3、5、7......个数值，列方向也取1、3、5、7......
个资料即可完成2:1:1采样),行、列方向都进行2:1:1抽样后的Cb、Cr矩阵大小
变成原来的1/4(M/2*S/2)。
现在的结果是Y矩阵大小为M*S，Cb,Cr为M/2*S/2，至此已经完成了颜色转换及采样，接下来做第二步--DCT变换。
二、DCT变换
在DCT变换之前得做一些准备工作，由于DCT变换一次只能做64个资料。因此，首先得把Y、Cb、Cr矩阵分成一个8*8的小块；其次，由于离散余弦变换公式所能接受的数值范围是-128至127之间，因此还得把Y、Cb、Cr矩阵中的每个资料减去128。
分块是一项比较烦琐的工作，尤其是对于Y矩阵，Cb、Cr矩阵好办，直接
按顺序分成8*8块就可以了，而Y块为了保证4个8*8Y块对应一个8*8Cb块和一个8*8Cr块，不致打乱顺序，得先把Y矩阵分成一个个16*16的块，这样得到一个16*16Y块对应一个8*8Cb块和一个8*8Cr块，然后再把每一个16*16块分成四个8*8小块，图标如下：
上图中每一个标有数字的小块表示一个8*8块，小块中的数字表示分块以
后，该8*8小块在矩阵中的位置顺序。
二维DCT变换公式为：
其中x,y代表图像数据矩阵中的某个资料值的坐标位置
f(x,y)指图像数据矩阵中该点的资料值
u,v指经过DCT变换后矩阵中的某数值点的坐标位置，在这里u=x,v=y
F(u,v)指经过DCT变换后该坐标点的资料值。
当u=0,v=0时，C(u)C(v)=1.414/2
当u>0,v>0时，C(u)C(v)=1,经过变换后的数据值F(u,v)称为频率系数(或
称DFT系数)。
由于VB汇编中都无法实现二维DCT计算公式，所以只有把二维的变换变
成先做一维，再做另一维的变换，一维DCT变换公式见附图一。
在我的程序中，设一个大小为8的数组SL(8),元素计数从1开始，而不是通常的0开始，先读取一个8*8块的第一行资料值，赋给SL(8),对SL(8)进行一维DCT变换后得到一个新的SL(8)数组，再把SL(8)数组覆盖到原来的8*8块中相应的地方去。做完第一行后再做第二行，一直做完8行，一个8*8块的一维DC即告完成，然后再做列方向的第二维DCT变换，变换公式一样，只是由SL(8)取8*8块的行资料变成取列数值。做完后覆盖回原值，即得到一个8*8块的DFT系数块，再重复这两个过程做第二个8*8块......一直到做完全部8*8块（Y，Cb，Cr）。这样就得到Y、Cb、Cr的DFT系数矩阵，要做到DCT系数矩阵还得把DFT系数矩阵以一个系数矩阵，Y矩阵用一个系数矩阵，Cb、Cr合用另一个系数矩阵，这两个系数矩阵见附录二，这两个系数矩阵可以直接得到第三步的量化结果。
三、量化
本来，第二步骤最后得到的DFT系数矩阵，乘以一个给定的系数矩阵得到DCT系数矩阵，DCT系数矩阵再除以一个量化矩阵得到量化后的DCT系数矩阵，而我们把这两步合成一步来做，即由系数矩阵和量化矩阵得到一个新的系数矩阵，亦即附录所列的矩阵，这样，把DFT系数数矩阵乘以新的系数矩阵就直接得到量化后的DCT矩阵。
量化过程实质上是把亮度资料和色度资料由室域转变成频域并滤除高频分量的过程，由于人眼对高频分量不敏感，所以可以滤除高频分量，经过量化以后的每一个8*8数据块中，第一个元素数据值为直流分量，称为DC，其余63个资料为交流分量，称为AC。
用程序实现量化过程，除了系数矩阵的推算比较麻烦之外，其余的都比较简
单，读取Y矩阵中第一个8%块，与量化系数矩阵中对应的相乘，得到的值覆盖回原矩阵，然后做第二个8*8块，一直到做完全部8*8块，然后做CbCr矩阵的量化，用另外一个系数矩阵。
经过量化后的DCT系数矩阵，除DC值一般不为零外，AC系数大多是在零点附近的浮点数。
由于霍夫曼编码的对象是整数，所以在做霍夫曼编码之前，还得对量化后的
DCT系数矩阵进行取整。
利用VB中的Int函数可以实现的YcbCr的DCT量化后系数炬阵进行取整。
经过取整以后，每一个8*8块中，有大量的AC系数的值为0。为了把尽可能多的其值为0的AC系数串在一起，以利于第四步的AC编码及提高压缩比，还必须把YcbCr矩阵中的每一个8*8块中的64个元素进行"之"字形排序（这样就可以做到把尽可能多的0串在一起），其过程示意图如下：假设SB（1）-SB（64）为一个8*8块中的64个元素，元素位置计数从1开始：
箭头方向表示"之"字形排序以后原8%中元素的新的位置顺序。亦即经过 "之"字形排序以下，新的8*8块中的元素值如下：等式左边表示元素位置，等式右边表示元素值为排序之前某位置中的元素值：如SB（15）＝SB（5），则左边表示排序后的8*8块中第15个元素的值等于排序之前第5个元素的值。
SB(1)=SB(1) SB(2)=SB(2)
SB(3)=SB(9) SB(4)=SB(17)
SB(5)=SB(10) SB(6)=SB(3)
SB(7)=SB(4) SB(8)=SB(11)
SB(9)=SB(18) SB(10)=SB(25)
SB(11)=SB(33) SB(12)=SB(26)
SB(13)=SB(19) SB(14)=SB(12)
SB(15)=SB(5) SB(12)=SB(6)
......
SB(61)=SB(48) SB(62)=SB(56)
SB(63)=SB(63) SB(64)=SB(64)
用程序实现：再新开一个大小为64的数组SC，把SB的值赋给SC，再把以上等式右边的SB换成SC即可。例如SB（15）＝CB（5）换成SB（15）＝SC（5）。换完以后即完成一个8*8块的之字形变换，这里为什么要新开一个数组SC呢？这是因为：例如把SB（9）的值赋给排序后的SB（3），则SB(3)的值就被SB(9)的值覆盖掉了，到后面需把SB (3)的值赋给SB(6)，如果不新开一个数组保留元素的原来的值，则赋给SB（6）的就成了是SB（9）的值而不是SB（3）的值，故要新开一个数组的保留8*8块中的原值。
完成"之"字形排序之后，就开始做第四步工作：霍夫曼编码。
四、霍夫曼编码
前面说到变换后的一个8*8频率系数矩阵由一个DC值和63个AC值构成，编码时对DC值和AC值用不同的霍夫曼编码表，对亮度和色度也需用不同的霍夫曼编码表，所以必须使用四张不同的霍夫曼编码表，才能完成JPEG编码。
DC编码：在做DC编码之前，还必须对DC值进行脉冲差值运算，具体做法是在Y、Cb、Cr频率系数矩阵中，后一个8*8块的DC值减去前一个8*8块的DC作为后一个8*8块新的DC值，并保留后一个8*8块的DC原值，用于后一个8*8块的差值DC运算，亦即每次后一个8*8块的DC值减去的是第一个8*8块的原来DC值，而不是经运算后的差值。
DC编码＝霍夫曼识别码（或称标志码）＋DC差值二进制代码
下面给出Y、CbCr矩阵的DC差值霍夫曼编码表。
举例：例如Y矩阵中一个8*8块的DC值=10，二进制值为1010，码长为4，查第一张表得霍夫曼识别长为3，识别码为101，则其霍夫曼编码为1011010，
负数的编码是取其绝对值的反码进行编码，如DC＝-10，则其绝对值为10，
二进制仍为1010，反码为0101，霍夫曼编码为1010101。
AC编码：AC编码的原理和方法跟DC相似，所不同的是AC编码中多了一项RLE压缩编码，前面说到经过量化取整以后，有许多AC值为0，并经过"之"字形排序，把原可能多的0串行在一起。在这里RLE压缩编码的就是用一个数值表示为0的AC值前有几个AC值为0。例如，在[M，N]这一组RLE编码中，N表示不为0的AC值，M则表示在这不为0的AC值N之前0的个数，M最多只能为15，如果AC资料值N之前有17个AC值为0，则先以[15，0]代表有16个值为0，再以[1，N]表示N前有一个值为0，如果在某个AC资料值之后（该值不为0），所有AC值皆为0，则这串资料可以用[0，0]表示。
做完RLE压缩后，再对不为0的AC值进行霍夫曼编码，跟DC值一样进行，只不过用的是另两张霍夫曼编码表，完整的AC编码如下：
（完整的AC编码码串包括三部分：（1）的位置记录"0"的个数；（2）的位置为霍夫曼识别码；（3）的位置的AC值的二进制代码值）这样的一个码串才算是一个完整的AC霍夫曼码串。
做到这里就完成了整个JPEG压缩过程。
在整个压缩过程中，最难的部分就是第四步，霍夫曼编码用程序实现起
非常繁琐，必须判断一个个DC（AC）的值，以及转换成二进制代码后的码长，
再去对照霍夫曼编码表进行编码，比如对一个DC值编码，首先得先判断该DC
的值在哪段范围内，在某一段范围内的数值，其二进制代码长相等。并要让程
序知道该DC的值到底为多少，然后才能进行编码，（这里面还涉及到一个数
制转换过程，即由十进制数转换成二进制数）。
压缩过程完成以后，接下去要做的工作便是码串存贮，存贮时需要注意
的-点是在抽样过程中得到4个Y块对应1个Cb块一个Cr块的对应关系，在存贮时也要按这种对应关系存贮，即4个Y块的码串后面紧接着存放与其对应的一个Cb块的码串及一个Cr块的码串，循环往复，存完所有串。这样的4个Y8*8块+1个Cb块+1个Cr简称为MCU，是JPEG格式的最小存贮处理单元。
附录1：一维DCT变换公式
C1 C2 C3 C4 C5 C6 C7为常数
C1=0.9808
C2=0.9239
C3=0.8315
C4=0.7071
C5=0.5556
C6=0.3827
C7=0.1951
S18=SL(1)+SL(8)
S27=SL(2)+SL(7)
S36=SL(3)+SL(6)
S45=SL(4)+SL(5)
S1845=S18+S45
S2736=S27+S36
D18=SL(1)+SL(8)
D27=SL(2)+SL(7)
D36=SL(3)+SL(6)
D45=SL(4)+SL(5)
D1845=S18-S45
D2736=S27-S36
SL(1)=0.5*(C4*(S1845+S2736))
SL(2)=0.5*(C1*D18+C3*D27+C5*D36+C7*D45)
SL(3)=0.5*(C2*D1845+C6*D2736)
SL(4)=0.5*(C3*D18-C7*D27-C1*D36-C5*D45)
SL(5)=0.5*(C4*(S1845-S2736))
SL(6)=0.5*(C5*D18-C1*D27+C7*D36+C3*D45)
SL(7)=0.5*(C6*D1845-C2*D2736)
SL(8)=0.5*(C7*D18-C5*D27+C3*D36-C1*D45)
附录2：DCT量化系数矩阵
