JPEG编码解码（Matlab）

搜索了网上的JPEG的matlab实现方式，发现只有寥寥几个，几乎都是只实现了一半，要么就是哈夫曼编码没有实现，要么就是只算出了哈夫曼的码长计算了一下效率，但是没有实际编码。要么就是太难我看不懂（汗。。。）比如github上几位大神的作品。在此只是简单的按照[1]中的过程实现了一个简单的JPEG，没有任何优化，只求简单。
在此简述一下JPEG编码解码的过程中需要关注的点。

其中仍有许多简化的部分：
（1）只计算了一张灰度图的编码解码，如果要是RGB通道的图需要额外处理RGB到YUV的变化，然后再分别对YUV进行编码，本文所讨论的即主要对其中的Y分量进行了编码。但是这部分的内容网上很容易搜索到。
（2）只对固定的图片进行了操作（512*512），因为没有添加jpeg文件头，所以很多变量都是我直接指定的。可以添加文件头来增加对多种文件格式的JPEG。
（3）存在一个假设，即每一个block在Z字形编排后后边全是0，如果存在一组最后不全是0，最后一位如果是个数的话这个代码可能就不行了，需要特殊处理一下。(这个可能性很小，但是应该依然存在，我猜。。。)

问题1：DCT变换
离散余弦变换的实现有三种方式[1]，第一种是用的矩阵的形式，这个在[4]中也采用的这种方式。[3]中详细的介绍了DCT的原理，非常好评！！！
为了方便没有书的同学，在此简述一下方式一：

当进行离散余弦变换时，$ 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内置的功能，大大简化了代码。当然出于练习的原因，我还是自己实现了一下这两行代码：

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()函数生效。

问题2：Z字形编排
[4]中使用了matlab自带的函数来处理这个问题，变得异常简单。

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