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 二进制类型数据相关操做