用VB語言編程實現JPEG數據壓縮

周俊敏

    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量化系數矩陣
 

原文轉自：http://www.anti-gravitydesign.com