JPEG編碼解碼（Matlab）

時間 2019-12-06

原文原文鏈接

搜索了網上的JPEG的matlab實現方式，發現只有寥寥幾個，幾乎都是隻實現了一半，要麼就是哈夫曼編碼沒有實現，要麼就是隻算出了哈夫曼的碼長計算了一下效率，可是沒有實際編碼。要麼就是太難我看不懂（汗。。。）好比github上幾位大神的做品。在此只是簡單的按照[1]中的過程實現了一個簡單的JPEG，沒有任何優化，只求簡單。
在此簡述一下JPEG編碼解碼的過程當中須要關注的點。html

其中仍有許多簡化的部分：
（1）只計算了一張灰度圖的編碼解碼，若是要是RGB通道的圖須要額外處理RGB到YUV的變化，而後再分別對YUV進行編碼，本文所討論的即主要對其中的Y份量進行了編碼。可是這部分的內容網上很容易搜索到。
（2）只對固定的圖片進行了操做（512*512），由於沒有添加jpeg文件頭，因此不少變量都是我直接指定的。能夠添加文件頭來增長對多種文件格式的JPEG。
（3）存在一個假設，即每個block在Z字形編排後後邊全是0，若是存在一組最後不全是0，最後一位若是是個數的話這個代碼可能就不行了，須要特殊處理一下。(這個可能性很小，可是應該依然存在，我猜。。。)git

問題1：DCT變換
離散餘弦變換的實現有三種方式[1]，第一種是用的矩陣的形式，這個在[4]中也採用的這種方式。[3]中詳細的介紹了DCT的原理，很是好評！！！
爲了方便沒有書的同窗，在此簡述一下方式一：github

當進行離散餘弦變換時，$ Y = AXA^T $ , 其中A即爲下邊生成的變換矩陣。X是輸入樣本矩陣，Y是變換後的係數矩陣。
進行逆變換時：$ X = A^{T}XA $ 。
其中A的公式以下：$ A_{ij} = C_i\cos{\frac{(2j+1)i\pi}{2N}} $
其中：$ C_i = \sqrt{\frac1N} \text{(i=0)}, C_i = \sqrt{\frac2N} \text{(i>0)}$算法

在matlab中只須要一行代碼就能夠實現：T = dctmtx(8)，固然咱們也能夠本身建立本身的變換矩陣：函數

N = 8;
    T = zeros(N,N);
    for i = 1:N
        for j = 1:N
            T(i, j) = sqrt(2/N) * cos(((2*j-1)*(i-1)*pi)/(2*N));
        end
    end
    T(1,:) = T(1,:) / sqrt(2);

dctmtx使用的是矩陣的方法，兩種方法結果類似。僅供參考。
固然也可使用matlab提供的dct2()和idct2()函數。這兩個個函數的核心實際上是第三種方法。優化

獲得了變換矩陣，進行DCT操做。
[4]中使用了一種很是簡便的方式，BY=blkproc(Y,[8 8],'P1*x*P2',T,T'); 進行DCT，能夠說很是方便。而後再使用BY2=blkproc(BY,[8 8],'round(x./P1)',a);進行量化。這兩個函數使用了matlab內置的功能，大大簡化了代碼。固然出於練習的緣由，我仍是本身實現了一下這兩行代碼：ui

quantization=zeros(X,Y);
        for i = 1:8:X
            for j = 1:8:Y
                mask = input_data(i:i+7,j:j+7);
                DCT = T * double(mask) * T';
                quantization(i:i+7,j:j+7) = round(DCT./Luminance_Quantization_Table);
            end
        end

順便把量化也作了。
在IDCT的時候：編碼

data=zeros(X, Y);
        
        for i = 1:8:X
            for j = 1:8:Y
                mask = decoded_matrix(i:i+7,j:j+7);
                mask =  mask .* Luminance_Quantization_Table;
                data(i:i+7,j:j+7) = T' * mask * T;
            end
        end
        data = uint8(data);

能夠看到恢復的時候有一些損失。注意：必定要把數據轉換成uint8，才能讓imshow()函數生效。.net

問題2：Z字形編排
[4]中使用了matlab自帶的函數來處理這個問題，變得異常簡單。code

% order
    order = [1 9  2  3  10 17 25 18 11 4  5  12 19 26 33  ...
        41 34 27 20 13 6  7  14 21 28 35 42 49 57 50  ...
        43 36 29 22 15 8  16 23 30 37 44 51 58 59 52  ...
        45 38 31 24 32 39 46 53 60 61 54 47 40 48 55  ...
        62 63 56 64];

注意的是order的順序和Z字形的序號並不相同，這是爲了讓這個表來適應matlab的特性。
使用:

y = im2col(quantization, [8 8], 'distinct');
        xb = size(y,2);
        y = y(order,:);

第一個函數是把整個圖劃分紅8*8=64的若干列，即每一列對應了一個8*8的塊，參數distinct是塊與塊不重疊。若使用默認sliding參數即爲窗口是滑動的。在8*8的塊排列的時候是按列排列的，窗口移動的過程當中也是豎向移動的。（matlab的默認操做都是以列爲基礎的，更好的使用列向量的一些特徵）

注意：在這裏最好把這個生成好的Y矩陣保存起來，這樣在解碼的時候能夠和解碼後的矩陣進行一下比較，容易發現問題。

恢復的時候，須要使用order來生成一個「反Z字形」序列：

rev = zeros(1,64);
        for k = 1:length(order)
            rev(k) = find(order==k);
        end
        
        X = 512;
        Y = 512;
        decoded = decoded(rev,:);
        decoded_matrix = col2im(decoded, [8 8], [X Y], 'distinct');

問題3：熵編碼（後兩個部分幾乎都是本身想的，可能效率很低）
量化表採用了直接做爲變量的形式存了起來。這個在[3]中的github代碼裏有C語言版的，複製過來稍加改動便可支持matlab。
而對於Huffman編碼表則是把[2]的碼錶直接以txt的形式存了起來。而後使用如下語句:

[ac_RS, ~, ac_code] = textread('AC_Hoffman_coding_table.txt', '%s%s%s');
    [dc_RS, ~, dc_code] = textread('DC_Hoffman_coding_table.txt', '%s%s%s');

把txt中的內容讀進來，讀取結果直接是三列，並且Length相等，方便下一步操做。
當進行轉換的時候直接使用相似一種查表的方式（DC編碼方式爲例）：

...
            if dc==0
                % 特殊狀況
                dc_encode = '00';
            else
                % 先對SSSS進行編碼
                SSSS = floor(log2(abs(dc)))+1;
                SSSS_index = strcmp(dc_RS ,string(SSSS));
                SSSS_encode = cell2mat(dc_code(SSSS_index));
                ...
                % 對DIFF進行編碼
                if dc > 0
                    DIFF_encode = dec2bin(dc);
                elseif dc < 0
                    dc_b = abs(dc);
                    dc_d = bitcmp(uint16(dc_b));
                    DIFF_encode = dec2bin(dc_d);
                    DIFF_encode = DIFF_encode(end-SSSS+1:end);
            end

注意，在進行小於0的數字進行編碼時，採起的應該是反碼（書上說的是補碼，沒看懂），這一點在[3]中有很是詳盡的說明。[3]寫的是真的好!在使用bitcmp()取反碼時，要注意的是取完反碼要捨去多餘的位，注意dec2bin的第二個參數指的是「at least n bits」，是個大坑。因此仍是本身截取一下吧。
在AC編碼時幾乎和DC相同，有幾個問題也須要注意，第一個是兩種特殊狀況的判斷。
特殊狀況1時若是後邊全是0了，就能夠直接結束了：

if mask(j:end)==0
                    ac_encode = [ac_encode '1010'];
                    break
                end

特殊狀況2是超過15個0了也要特殊處理：

if mask(j)==0
                    zero_tot = zero_tot + 1;
                    if zero_tot == 16
                        ac_encode = [ac_encode '11111111001'];
                        zero_tot = 0;
                    end
                    continue
                end

而後就是要把RRRR和SSSS拼接起來，在進行查找：

S = floor(log2(abs(mask(j))))+1;
                R_S = [dec2hex(zero_tot) '/' dec2hex(S)];

問題4：解碼
解碼時先對熵編碼進行解碼，而後和剛纔存好的y進行比較，看看是否是同樣。而後再恢復Z字形，再IDCT（這兩個前邊已經寫過了）。
熵編碼解碼的時候沒次生成一個64行一列的mask，而後拼接到一塊兒。解碼的時候先拿着碼串取AC_cod和DC_code中匹配，能匹配成功再進行下一步

index = strcmp(dc_code, tmp);
                
                if any(index)
                    SSSS = cell2mat(dc_RS(index));
                    SSSS = str2double(SSSS);
                    
                    if SSSS == 0
                        DIFF = '0';
                    else
                        DIFF = input_data(p_end+1:p_end+SSSS);
                    end

使用any能夠避免一些坑。也要注意0的處理。
注意首位是0的時候，要進行反碼操做：

if DIFF(1) == '0'
                        dc_d = bin2dec(DIFF);
                        dc_b = bitcmp(uint16(dc_d));
                        dc = -double(mod(dc_b, 2^SSSS));
                    else
                        dc = bin2dec(DIFF);
                    end

dc = -double(mod(dc_b, 2^SSSS)); 這一句若是不加double會出現餘數爲1可是取負號就成了0的詭異狀況。。。難以解釋。
對AC解碼時須要注意的是對「0/0」和「F/0」分別處理一下就好。

Reference
[1] 《多媒體技術基礎（第4版）》林福宗清華大學出版社 131-140
[2] JPEG 標準推薦的亮度、色度DC、AC Huffman 編碼表
[3] JPEG算法解密
[4] 使用Matlab實現JPEG壓縮 (github)
[5] Matlab 二進制類型數據相關操做