JPEG原理分析及JPEG解码器的调试

JPEG原理分析及JPEG解码器的调试

一、JPEG格式介绍

$\text{§}$ 1.1 JPEG简介

JPEG 是 Joint Photographic Exports Group 的英文缩写，中文称之为联合图像专家小组。该小组隶属于 ISO 国际标准化组织，主要负责定制静态数字图像的编码方法，即所谓的 JPEG算法。
JPEG图像压缩算法能够在提供良好的压缩性能的同时，具有比较好的重建质量，被广泛应用于图像、视频处理领域，网站上80%的图像都采用了JPEG压缩标准。

$\text{§}$ 1.2 JPEG编码流程

过程解析：

1. 将传入的图像的亮度分量全部进行128电平偏置，即$pixe{l}_{val}=pixe{l}_{val}-128$，以此将每个像素点的取值范围从$[0,255]$变换到$[-128,127]$。这个部分其实是一个历史遗留问题，但是也可以发现这一操作使我们的算法系统设计更为便利。
2. 进行$8\times 8$的DCT变换，实现能量集中和去相关，降低空间冗余度。
3. 进行均匀量化，这一过程需要使用由人眼特性制作的量化表，特点是低频细量化，高频粗量化。
4. 对直流系数进行差分和VLC编码，对交流系数使用zig-zag扫描和游程编码，再进行VLC编码。

$\text{§}$ 1.2.1 颜色空间转换

JPEG采用的是YCbCr颜色空间，所以在传入RGB图片前，需要先进行色彩空间转换,YCrCb颜色空间中，Y 代表亮度，Cr,Cb则代表色度和饱和度(也有人将 Cb,Cr 两者统称为色度)，三者通常以 Y,U,V 来表示，即用 U 代表 Cb，用 V 代表 Cr。
$$\begin{gathered} Y=0.299R+0.587G+0.114B \\ Cb=-0.1687R-0.3313G+0.5B+128 \\ Cr=0.5R=0.418G-0.0813B+128 \end{gathered}$$

$\text{§}$ 1.2.2 采样

人眼对于亮度信号的敏感度比对色度信号的敏感度更高，因此一般我们选用的图像都是4:2:0或者4:2:2的采样格式。

$\text{§}$ 1.2.3 电平偏置

如上述过程解析中所说，需要对传入的每个像素值减去128。

$\text{§}$ 1.2.4 $8\times 8$DCT变换

DCT 变换是对$8\times 8$ 的子块进行处理的，因此，在进行 DCT 变换之前必须把源图象数据进行分块。源图像中每点的 3个分量是交替出现的，先要把这 3 个分量分开，存放到 3 张表中去。然后由左及右，由上到下依次读取$8\times 8$的子块，存放在长度为 64 的表中，即可以进行 DCT 变换。这里注意，如果图像的宽高不是8的倍数，需要进行填充，通常采取填充相邻像素值的方法进行，最大限度避免出现高频分量。

DCT变换（离散余弦变换）使傅里叶变换相关的一种变换，通过DCT变换可以使图像能量集中，同时去除相关性，降低冗余度。

$\text{§}$ 1.2.5 量化

对进行DCT变换后的系数进行量化，会有一张量化表规定每个DCT系数的量化步长。

$\text{§}$ 1.2.6 DC系数编码

直流系数$f(0,0)$反映了该图像包含的全部直流成分，一般值会很大，又因为通常两个相邻子图像块之间具有较强的相关性，所以可以引入DPCM差分预测编码。即：对本像素块的直流系数和前一像素块的直流系数的差值进行编码。

$\text{§}$ 1.2.7 AC系数编码

经过DCT变换后，大部分的非零系数都集中在了DCT矩阵的左上角，通过Z字型读取后我们可以尽可能将0系数集中在一起，方便后续的游程编码。这一部分的实现在代码中体现为类似八皇后问题的思路：

zig-zag扫描后，进行游程编码（RLC）。在JPEG中，游程编码的格式规定为：$(run,level)$。

• 表示连续run个0，后面跟一个值为level的系数。
• run最多15个，用4位表示XXXX。
• level类似DC，分成16个类别，用4位表示类别号RRRR，并支持类内索引。
• 对(XXXX,RRRR)进行Huffman编码。
• 对类内索引进行定长码编码

$\text{§}$ 1.3 JPEG文件格式

$\text{§}$ 1.3.1 JPEG语法结构

JPEG在文件中以Segment的形式组织，具有以下特点：

• 均以0xFF开始，后跟1字节的标记标识符和2字节的标记长度以及该标记所对应的payloada。
• 标记长度部分高位在前，低位在后，不包含该标记的头两个字节。
• 熵编码部分的数据在0xFF后由编码器插入0x00，解码器解码时跳过此字节不予处理。
• SOI和EOI标记没有payload。

参数说明：

• SOI：Start Of Image表示图像的开始，固定值为FFD8。
• SOF0：帧图像的开始，记录每一帧图像的数据长度、样本数据的位数、图像的宽高、颜色分量数、颜色分量信息（分量ID（Y、U、v）、采样因子（4：4：4、4：2：2、4：2：0）、量化表ID）。
• DQT：DCT量化表，记录量化表长度、量化精度、量化表ID、表项（长度为64bit(8位精度)，记录了8*8DCT变换后每个像素的量化步长，由于DC、AC、亮度、色度使用不同的量化编，所有量化表最多有4个）。
• DHT0：定义Huffman码表。表长度、表ID（0：亮度 1：色度）、表类型（0：直流 1：交流）不同位数的码字数量（16字节分别记录了长度为1到16的码字的个数）、权值。
• APP0：应用程序保留标记 （版本参数信息）。
• DC表：权值的大小直流分量数值的二进制位数，读取后经过查表查得对应的DC值。权值的字节数为DC经DPCM编码后码字个数的总和。
• AC表：权值的高四位表示当前数值前面有多少个0，低4为表示交流分量数值的二进制位数。
• SOS：扫描开始标记，记录了数据长度、颜色分量数（与SOF0相同）、颜色分量信息（颜色分量ID：１２３对应ＹＵＶ）、表号：（高位为直流系数使用的Huffman表数、低位为交流系数使用的Huffman表数）、压缩图像数据。
• EOI：End Of Image，图像结束标记。

注意：

1. JPEG由若干个必不可少的标记顺序连接构成整个文件。
2. 文件一定以0xFFD8开始，即表示图像开始的SOI标记。
3. 文件一定以0xFFD9结束，表示图像结束的EOI标记。

$\text{§}$ 1.3.2 APP0 Segment

有九个字段

字段	字节数	含义
数据长度	2字节	1-9共9个字段的总长度
标识符	5字节	固定值0x4A46494600，即字符串“JFIF0”
版本号	2字节	一般是0x0102，表示JFIF的版本号1.2
X和Y的密度单位	1字节	只有三个值可选 0：无单位；1：点数/英寸；2：点数/厘米
X方向像素密度	2字节	X方向像素密度
Y方向像素密度	2字节	Y方向像素密度
缩略图水平像素数	1字节	缩略图水平像素数目
缩略图垂直像素数	1字节	缩略图垂直像素数目
缩略图RGB位图	3字节的倍数	缩略RGB位图，n为像素数

$\text{§}$ 1.3.3 量化表DQT

• 一般是两个量化表，即亮度和色度各一张。
• 以0xFFDB开始：
  • 量化表长度，一般是0043或者0084。
  • 量化表信息（1字节）
    • bit 0-3 ：QT号，取值只能是0-3。
    • bit 4-7：QT精度，0为8比特，否则是16比特。
  • 量化表的实际数据，按照Z字形保存量化表内$8\times 8$的数据。

$\text{§}$ 1.3.4 帧图像开始SOF0

以0xFFC0开始

字段	字节数	含义
SOF长度	2字节	固定为0x0011
精度	1字节	每个颜色分量每个像素的位数，通常为8
图像高度	2字节	以像素数表示图像的高度
图像宽度	2字节	以像素数表示图像的宽度
颜色分量数	1字节	通常是3
颜色分量信息	颜色分量数x3字节，一般是9字节	颜色分量ID 1字节 水平/垂直采样因子 1字节 高四位为水平采样因子，低四位为垂直采样因子 量化表 1字节 当前分量所使用量化表的ID

$\text{§}$ 1.3.5 Huffman表

数值为0xFFC4。

• Huffman表的长度
• 类型（AC或者DC）
• 索引index
• 位表bit table
• 值表value table

其中Huffman表的类型只有四种：

• 0x00表示DC直流0号表
• 0x01表示DC直流1号表
• 0x10表示AC交流0号表
• 0x11表示AC交流1号表

举例说明：如下是一个Huffman码表

• 黑色部分包含表头和长度（0x001D）
• 红色部分是表的ID和类型，这里的00表示此部分数据描述的是DC直流0号表。
• 蓝色部分（16字节），为不同位数的码字的数目。这16个数值的实际意义是：没有1位的码字；有3个2位的码字；有1个3位的码字…以此类推。
• 绿色部分（10字节）：为编码内容，通过蓝色数据部分我们知道Huffman树中一共有10个叶子节点，所以这一部分的数据长度就是10字节。表示每个叶子节点从小到大(所谓从小到大实际上是指的码长和出现顺序)排序后，对应的权值(权值对直流和交流系数的含义不同，这一点在后面会说明)。

建立Huffman表

建立Huffman表的方式其实很简单：

具体方法为：

• 第一个码字一定是0
  • 如果第一个码字位为1，则码字为0；
  • 如果第一个码字位为2，则为00；
  • 以此类推…
• 从第二个码字开始
  • 如果它和它前面的码字位数相同，则当前码字为它前面的码字加1；
  • 如果它的位数比它前面的码字位数大，则当前码字是前面的码字加1后再在后边添若干个0，直至满足位数长度为止。

我们按照上面给出的例子建立一个Huffman码表：

序号	码长	码字	权重
1	2	00	4
2	2	01	5
3	2	10	6
4	3	110	3
5	4	1110	2
6	5	11110	1
7	6	111110	0
8	7	1111110	9
9	8	11111110	7
10	9	111111110	8

下面解释一下在直流和交流中权重的作用：

直流

对于直流来说，权值是解码时需要额外读入的bit位数。这个再次读入的数通过查表来得到真正的码值。
比如比特流：0110101011
读入01后码字结束，查表可以得到额外比特数为5，所以需要额外读5位数，此处为10101，译码可以得到21，所以直流系数为21。要注意的是，直流系数是差分编码之后得到的。

实际中的过程是把所有的颜色分量单元按颜色分量（Y、Cr、Cb）分类。每一种颜色分量内，相邻的两个颜色分量单元的直流变量是以差分来编码的。也就是说，通过Huffman码表解码出来的直流变量数值只是当前颜色分量单元的实际直流变量减去前一个颜色分量单元的实际直流变量。也就是说，当前直流变量要通过前一个颜色分量单元的实际（非解码）直流分量来校正，即$D{C}_n=D{C}_{n-1}+Diff$，这里的$Diff$就是我们解码得到的直流系数，如果此时的颜色分量单元是第一个单元，则解码出来的直流系数就是真正的直流变量值。

交流

对于交流系数，用交流哈夫曼树/表查得该码字对应的权值。权值的高4位表示当前数值前面有多少个连续的零，低4位表示该交流分量数值的二进制位数，也就是接下来需要读入的位数。

例如，权值0x31可以表示为(3,1)。表明交流系数前面有3个0，此外交流系数的具体值还需要再读入1bit的码字才能得到。

二、JPEG代码分析

$\text{§}$ 2.1 代码架构

$\text{§}$ 2.1.1 main函数

接受输入输出文件名称参数，打开TRACEFILE。通过输入参数选择需要输出的文件格式，此实验为YUV420，同时调用convert_one_image()。

int main(int argc, char* argv[])
{
    int output_format = TINYJPEG_FMT_YUV420P;
    char* output_filename, * input_filename;
    clock_t start_time, finish_time;
    unsigned int duration;
    int current_argument;
    int benchmark_mode = 0;
#if TRACE
    p_trace = fopen(TRACEFILE, "w");
    if (p_trace == NULL)
    {
        printf("trace file open error!");
    }
#endif
    if (argc < 3)
        usage();

    current_argument = 1;
    while (1)
    {
        if (strcmp(argv[current_argument], "--benchmark") == 0)
            benchmark_mode = 1;
        else
            break;
        current_argument++;
    }

    if (argc < current_argument + 2)
        usage();

    input_filename = argv[current_argument];
    if (strcmp(argv[current_argument + 1], "yuv420p") == 0)
        output_format = TINYJPEG_FMT_YUV420P;
    else if (strcmp(argv[current_argument + 1], "rgb24") == 0)
        output_format = TINYJPEG_FMT_RGB24;
    else if (strcmp(argv[current_argument + 1], "bgr24") == 0)
        output_format = TINYJPEG_FMT_BGR24;
    else if (strcmp(argv[current_argument + 1], "grey") == 0)
        output_format = TINYJPEG_FMT_GREY;
    else
        exitmessage("Bad format: need to be one of yuv420p, rgb24, bgr24, grey\n");
    output_filename = argv[current_argument + 2];

    start_time = clock();

    if (benchmark_mode)
        load_multiple_times(input_filename, output_filename, output_format);
    else
        convert_one_image(input_filename, output_filename, output_format);

    finish_time = clock();
    duration = finish_time - start_time;
    snprintf(error_string, sizeof(error_string), "Decoding finished in %u ticks\n", duration);
#if TRACE
    fclose(p_trace);
#endif
    return 0;
}

$\text{§}$ 2.1.2 convert_one_image

打开输入输出文件，初始化jdec结构体，获得文件参数信息，主要调用tinyjpeg_parse_header()解码jpeg图像，最后调用write_yuv()写入yuv文件。

/**
 * Load one jpeg image, and decompress it, and save the result.
 * 加载一个JPEG图像，并进行解压缩和存储
 */
int convert_one_image(const char *infilename, const char *outfilename, int output_format)
{
  FILE *fp;
  unsigned int length_of_file;
  unsigned int width, height;
  unsigned char *buf;
  struct jdec_private *jdec;
  unsigned char *components[3];

  /* Load the Jpeg into memory 加载图像*/
  fp = fopen(infilename, "rb");
  if (fp == NULL)
    exitmessage("Cannot open filename\n");
  length_of_file = filesize(fp);
  buf = (unsigned char *)malloc(length_of_file + 4);
  if (buf == NULL)
    exitmessage("Not enough memory for loading file\n");
  fread(buf, length_of_file, 1, fp);
  fclose(fp);

  /* Decompress it 初始化解码器*/
  jdec = tinyjpeg_init();
  if (jdec == NULL)
    exitmessage("Not enough memory to alloc the structure need for decompressing\n");
  // 解码头部信息
  if (tinyjpeg_parse_header(jdec, buf, length_of_file)<0)
    exitmessage(tinyjpeg_get_errorstring(jdec));

  /* Get the size of the image */
  tinyjpeg_get_size(jdec, &width, &height);

  snprintf(error_string, sizeof(error_string),"Decoding JPEG image...\n");
  // 开始解码
  if (tinyjpeg_decode(jdec, output_format) < 0)
    exitmessage(tinyjpeg_get_errorstring(jdec));

  /* 
   * Get address for each plane (not only max 3 planes is supported), and
   * depending of the output mode, only some components will be filled 
   * RGB: 1 plane, YUV420P: 3 planes, GREY: 1 plane
   */
  tinyjpeg_get_components(jdec, components);

  /* Save it */
  switch (output_format)
   {
    case TINYJPEG_FMT_RGB24:
    case TINYJPEG_FMT_BGR24:
      write_tga(outfilename, output_format, width, height, components);
      break;
    case TINYJPEG_FMT_YUV420P:
      write_yuv(outfilename, width, height, components);
      break;
    case TINYJPEG_FMT_GREY:
      write_pgm(outfilename, width, height, components);
      break;
   }

  /* Only called this if the buffers were allocated by tinyjpeg_decode() */
  tinyjpeg_free(jdec);
  /* else called just free(jdec); */

  free(buf);
  return 0;
}

$\text{§}$ 2.1.3 tinyjpeg_decode

/**
 * Decode and convert the jpeg image into @pixfmt@ image
 * 
 * Note: components will be automaticaly allocated if no memory is attached.
 */
int tinyjpeg_decode(struct jdec_private *priv, int pixfmt)
{
  unsigned int x, y, xstride_by_mcu, ystride_by_mcu;
  unsigned int bytes_per_blocklines[3], bytes_per_mcu[3];
  decode_MCU_fct decode_MCU;
  const decode_MCU_fct *decode_mcu_table;
  const convert_colorspace_fct *colorspace_array_conv;
  convert_colorspace_fct convert_to_pixfmt;

  if (setjmp(priv->jump_state))
    return -1;

  /* To keep gcc happy initialize some array */
  bytes_per_mcu[1] = 0;
  bytes_per_mcu[2] = 0;
  bytes_per_blocklines[1] = 0;
  bytes_per_blocklines[2] = 0;

  decode_mcu_table = decode_mcu_3comp_table;
  switch (pixfmt) {
     case TINYJPEG_FMT_YUV420P:
       colorspace_array_conv = convert_colorspace_yuv420p;
       if (priv->components[0] == NULL)
	 priv->components[0] = (uint8_t *)malloc(priv->width * priv->height);
       if (priv->components[1] == NULL)
	 priv->components[1] = (uint8_t *)malloc(priv->width * priv->height/4);
       if (priv->components[2] == NULL)
	 priv->components[2] = (uint8_t *)malloc(priv->width * priv->height/4);
       bytes_per_blocklines[0] = priv->width;
       bytes_per_blocklines[1] = priv->width/4;
       bytes_per_blocklines[2] = priv->width/4;
       bytes_per_mcu[0] = 8;
       bytes_per_mcu[1] = 4;
       bytes_per_mcu[2] = 4;
       break;

     case TINYJPEG_FMT_RGB24:
       colorspace_array_conv = convert_colorspace_rgb24;
       if (priv->components[0] == NULL)
	 priv->components[0] = (uint8_t *)malloc(priv->width * priv->height * 3);
       bytes_per_blocklines[0] = priv->width * 3;
       bytes_per_mcu[0] = 3*8;
       break;

     case TINYJPEG_FMT_BGR24:
       colorspace_array_conv = convert_colorspace_bgr24;
       if (priv->components[0] == NULL)
	 priv->components[0] = (uint8_t *)malloc(priv->width * priv->height * 3);
       bytes_per_blocklines[0] = priv->width * 3;
       bytes_per_mcu[0] = 3*8;
       break;

     case TINYJPEG_FMT_GREY:
       decode_mcu_table = decode_mcu_1comp_table;
       colorspace_array_conv = convert_colorspace_grey;
       if (priv->components[0] == NULL)
	 priv->components[0] = (uint8_t *)malloc(priv->width * priv->height);
       bytes_per_blocklines[0] = priv->width;
       bytes_per_mcu[0] = 8;
       break;

     default:
#if TRACE
		 fprintf(p_trace,"Bad pixel format\n");
		 fflush(p_trace);
#endif
       return -1;
  }

  xstride_by_mcu = ystride_by_mcu = 8;
  if ((priv->component_infos[cY].Hfactor | priv->component_infos[cY].Vfactor) == 1) {
     decode_MCU = decode_mcu_table[0];
     convert_to_pixfmt = colorspace_array_conv[0];
#if TRACE
     fprintf(p_trace,"Use decode 1x1 sampling\n");
	 fflush(p_trace);
#endif
  } else if (priv->component_infos[cY].Hfactor == 1) {
     decode_MCU = decode_mcu_table[1];
     convert_to_pixfmt = colorspace_array_conv[1];
     ystride_by_mcu = 16;
#if TRACE
     fprintf(p_trace,"Use decode 1x2 sampling (not supported)\n");
	 fflush(p_trace);
#endif
  } else if (priv->component_infos[cY].Vfactor == 2) {
     decode_MCU = decode_mcu_table[3];
     convert_to_pixfmt = colorspace_array_conv[3];
     xstride_by_mcu = 16;
     ystride_by_mcu = 16;
#if TRACE 
	 fprintf(p_trace,"Use decode 2x2 sampling\n");
	 fflush(p_trace);
#endif
  } else {
     decode_MCU = decode_mcu_table[2];
     convert_to_pixfmt = colorspace_array_conv[2];
     xstride_by_mcu = 16;
#if TRACE	
     fprintf(p_trace,"Use decode 2x1 sampling\n");
	 fflush(p_trace);
#endif
  }

  resync(priv);

  /* Don't forget to that block can be either 8 or 16 lines */
  bytes_per_blocklines[0] *= ystride_by_mcu;
  bytes_per_blocklines[1] *= ystride_by_mcu;
  bytes_per_blocklines[2] *= ystride_by_mcu;

  bytes_per_mcu[0] *= xstride_by_mcu/8;
  bytes_per_mcu[1] *= xstride_by_mcu/8;
  bytes_per_mcu[2] *= xstride_by_mcu/8;

  /* Just the decode the image by macroblock (size is 8x8, 8x16, or 16x16) */
  // 分块进行编码
  for (y=0; y < priv->height/ystride_by_mcu; y++)
   {
     //trace("Decoding row %d\n", y);
     priv->plane[0] = priv->components[0] + (y * bytes_per_blocklines[0]);
     priv->plane[1] = priv->components[1] + (y * bytes_per_blocklines[1]);
     priv->plane[2] = priv->components[2] + (y * bytes_per_blocklines[2]);
     for (x=0; x < priv->width; x+=xstride_by_mcu)
      {
	decode_MCU(priv);
	convert_to_pixfmt(priv);
	priv->plane[0] += bytes_per_mcu[0];
	priv->plane[1] += bytes_per_mcu[1];
	priv->plane[2] += bytes_per_mcu[2];
	if (priv->restarts_to_go>0)
	 {
	   priv->restarts_to_go--;
	   if (priv->restarts_to_go == 0)
	    {
	      priv->stream -= (priv->nbits_in_reservoir/8);
	      resync(priv);
	      if (find_next_rst_marker(priv) < 0)
		return -1;
	    }
	 }
      }
   }
#if TRACE
  fprintf(p_trace,"Input file size: %d\n", priv->stream_length+2);
  fprintf(p_trace,"Input bytes actually read: %d\n", priv->stream - priv->stream_begin + 2);
  fflush(p_trace);
#endif

  return 0;
}

$\text{§}$ 2.1.4 tinyjpeg_parse_header

控制指针移动，调用parse_JFIF()

/**
 * Initialize the tinyjpeg object and prepare the decoding of the stream.
 * 初始化并准备后续解码
 * Check if the jpeg can be decoded with this jpeg decoder.
 * 检查我们给出的JPEG图像是否可以被JPEG解码器解码
 * Fill some table used for preprocessing.
 */
int tinyjpeg_parse_header(struct jdec_private *priv, const unsigned char *buf, unsigned int size)
{
  int ret;

  /* Identify the file */
  if ((buf[0] != 0xFF) || (buf[1] != SOI))
    snprintf(error_string, sizeof(error_string),"Not a JPG file ?\n");

  priv->stream_begin = buf+2;
  priv->stream_length = size-2;
  priv->stream_end = priv->stream_begin + priv->stream_length;

  ret = parse_JFIF(priv, priv->stream_begin);

  return ret;
}

$\text{§}$ 2.1.5 parse_JFIF

循环调用parse_SOF、parse_DQT、parse_SOS、parse_DHT、parse_DRI函数，直到发现SOS块为止。

static int parse_JFIF(struct jdec_private *priv, const unsigned char *stream)
{
  int chuck_len;
  int marker;
  int sos_marker_found = 0;
  int dht_marker_found = 0;
  const unsigned char *next_chunck;

  /* Parse marker */
  // 循环调用，直到找到sos
  while (!sos_marker_found)
   {
     if (*stream++ != 0xff)
       goto bogus_jpeg_format;
     /* Skip any padding ff byte (this is normal) */
     while (*stream == 0xff)
       stream++;

     marker = *stream++;
     chuck_len = be16_to_cpu(stream);
     next_chunck = stream + chuck_len;
     switch (marker)
      {
       case SOF:
	 if (parse_SOF(priv, stream) < 0)
	   return -1;
	 break;
       case DQT:
	 if (parse_DQT(priv, stream) < 0)
	   return -1;
	 break;
       case SOS:
	 if (parse_SOS(priv, stream) < 0)
	   return -1;
	 sos_marker_found = 1;
	 break;
       case DHT:
	 if (parse_DHT(priv, stream) < 0)
	   return -1;
	 dht_marker_found = 1;
	 break;
       case DRI:
	 if (parse_DRI(priv, stream) < 0)
	   return -1;
	 break;
       default:
#if TRACE
	fprintf(p_trace,"> Unknown marker %2.2x\n", marker);
	fflush(p_trace);
#endif
	 break;
      }

     stream = next_chunck;
   }

  if (!dht_marker_found) {
#if TRACE
	  fprintf(p_trace,"No Huffman table loaded, using the default one\n");
	  fflush(p_trace);
#endif
    build_default_huffman_tables(priv);
  }

#ifdef SANITY_CHECK
  if (   (priv->component_infos[cY].Hfactor < priv->component_infos[cCb].Hfactor)
      || (priv->component_infos[cY].Hfactor < priv->component_infos[cCr].Hfactor))
    snprintf(error_string, sizeof(error_string),"Horizontal sampling factor for Y should be greater than horitontal sampling factor for Cb or Cr\n");
  if (   (priv->component_infos[cY].Vfactor < priv->component_infos[cCb].Vfactor)
      || (priv->component_infos[cY].Vfactor < priv->component_infos[cCr].Vfactor))
    snprintf(error_string, sizeof(error_string),"Vertical sampling factor for Y should be greater than vertical sampling factor for Cb or Cr\n");
  if (   (priv->component_infos[cCb].Hfactor!=1) 
      || (priv->component_infos[cCr].Hfactor!=1)
      || (priv->component_infos[cCb].Vfactor!=1)
      || (priv->component_infos[cCr].Vfactor!=1))
    snprintf(error_string, sizeof(error_string),"Sampling other than 1x1 for Cr and Cb is not supported");
#endif

  return 0;
bogus_jpeg_format:
#if TRACE
  fprintf(p_trace,"Bogus jpeg format\n");
  fflush(p_trace);
#endif
  return -1;
}

$\text{§}$ 2.1.6 parse_DQT

调用build_quantization_table函数创建量化表，解码获得量化表信息。

static int parse_DQT(struct jdec_private *priv, const unsigned char *stream)
{
  int qi;
  float *table;
  const unsigned char *dqt_block_end;
#if TRACE
  fprintf(p_trace,"> DQT marker\n");
  fflush(p_trace);
#endif
  dqt_block_end = stream + be16_to_cpu(stream);
  stream += 2;	/* Skip length */

  while (stream < dqt_block_end)
   {
     qi = *stream++;
#if SANITY_CHECK
     if (qi>>4)
       snprintf(error_string, sizeof(error_string),"16 bits quantization table is not supported\n");
     if (qi>4)
       snprintf(error_string, sizeof(error_string),"No more 4 quantization table is supported (got %d)\n", qi);
#endif
     table = priv->Q_tables[qi];
     build_quantization_table(table, stream);
     stream += 64;
   }
#if TRACE
  fprintf(p_trace,"< DQT marker\n");
  fflush(p_trace);
#endif
  return 0;
}

$\text{§}$ 2.1.7 parse_DHT

解码获得Huffman码表信息，通过build_huffman_table创建Huffman表。

static int parse_DHT(struct jdec_private *priv, const unsigned char *stream)
{
  unsigned int count, i;
  unsigned char huff_bits[17];
  int length, index;

  length = be16_to_cpu(stream) - 2;
  stream += 2;	/* Skip length */
#if TRACE
  fprintf(p_trace,"> DHT marker (length=%d)\n", length);
  fflush(p_trace);
#endif

  while (length>0) {
     index = *stream++;

     /* We need to calculate the number of bytes 'vals' will takes */
     huff_bits[0] = 0;
     count = 0;
     for (i=1; i<17; i++) {
	huff_bits[i] = *stream++;
	count += huff_bits[i];
     }
#if SANITY_CHECK
     if (count >= HUFFMAN_BITS_SIZE)
       snprintf(error_string, sizeof(error_string),"No more than %d bytes is allowed to describe a huffman table", HUFFMAN_BITS_SIZE);
     if ( (index &0xf) >= HUFFMAN_TABLES)
       snprintf(error_string, sizeof(error_string),"No more than %d Huffman tables is supported (got %d)\n", HUFFMAN_TABLES, index&0xf);
#if TRACE
     fprintf(p_trace,"Huffman table %s[%d] length=%d\n", (index&0xf0)?"AC":"DC", index&0xf, count);
	 fflush(p_trace);
#endif
#endif

     if (index & 0xf0 )
       build_huffman_table(huff_bits, stream, &priv->HTAC[index&0xf]);
     else
       build_huffman_table(huff_bits, stream, &priv->HTDC[index&0xf]);

     length -= 1;
     length -= 16;
     length -= count;
     stream += count;
  }
#if TRACE
  fprintf(p_trace,"< DHT marker\n");
  fflush(p_trace);
#endif
  return 0;
}

$\text{§}$ 2.2 结构体解析

$\text{§}$ 2.2.1 huffman_table

用于存储Huffman码表，包含快速查找表lookup，存储码字长度的表code_size，慢速表slowtable。

struct huffman_table
{
  /* Fast look up table, using HUFFMAN_HASH_NBITS bits we can have directly the symbol,
   * if the symbol is <0, then we need to look into the tree table */
  short int lookup[HUFFMAN_HASH_SIZE];
  /* code size: give the number of bits of a symbol is encoded */
  unsigned char code_size[HUFFMAN_HASH_SIZE];
  /* some place to store value that is not encoded in the lookup table 
   * FIXME: Calculate if 256 value is enough to store all values
   */
  uint16_t slowtable[16-HUFFMAN_HASH_NBITS][256];
};

$\text{§}$ 2.2.2 component

存储DCT变换后的数值及指向使用的量化表与霍夫曼码表

struct component 
{
  unsigned int Hfactor;
  unsigned int Vfactor;
  float *Q_table;		/* Pointer to the quantisation table to use */
  struct huffman_table *AC_table;
  struct huffman_table *DC_table;
  short int previous_DC;	/* Previous DC coefficient */
  short int DCT[64];		/* DCT coef */
#if SANITY_CHECK
  unsigned int cid;
#endif
};

$\text{§}$ 2.2.3 jdec_private

是一个完整的JPEG码流块的定义，其中也包含了components和huffman_table的结构体。

struct jdec_private
{
  /* Public variables */
  uint8_t *components[COMPONENTS];
  unsigned int width, height;	/* Size of the image */
  unsigned int flags;

  /* Private variables */
  const unsigned char *stream_begin, *stream_end;
  unsigned int stream_length;

  const unsigned char *stream;	/* Pointer to the current stream */
  unsigned int reservoir, nbits_in_reservoir;

  struct component component_infos[COMPONENTS];
  float Q_tables[COMPONENTS][64];		/* quantization tables */
  struct huffman_table HTDC[HUFFMAN_TABLES];	/* DC huffman tables   */
  struct huffman_table HTAC[HUFFMAN_TABLES];	/* AC huffman tables   */
  int default_huffman_table_initialized;
  int restart_interval;
  int restarts_to_go;				/* MCUs left in this restart interval */
  int last_rst_marker_seen;			/* Rst marker is incremented each time */

  /* Temp space used after the IDCT to store each components */
  uint8_t Y[64*4], Cr[64], Cb[64];

  jmp_buf jump_state;
  /* Internal Pointer use for colorspace conversion, do not modify it !!! */
  uint8_t *plane[COMPONENTS];

};

三、实验结果

$\text{§}$ 3.1 JPEG to YUV

实验目的：将输入的JPEG文件转换为可以在YUVReview中查看的YUV文件。
为了方便实验，我们选取的图片最好是宽高为8的倍数，这里我们选取了分辨率为$1024\times 1024$的图片：

为了输出为YUV格式的文件，我们找到需要调用的函数：write_yuv，在原有代码的基础上添加组合为yuv文件的代码：

static void write_yuv(const char *filename, int width, int height, unsigned char **components)
{
  FILE *F;
  char temp[1024];

  snprintf(temp, 1024, "%s.Y", filename);
  F = fopen(temp, "wb");
  fwrite(components[0], width, height, F);
  fclose(F);
  snprintf(temp, 1024, "%s.U", filename);
  F = fopen(temp, "wb");
  fwrite(components[1], width*height/4, 1, F);
  fclose(F);
  snprintf(temp, 1024, "%s.V", filename);
  F = fopen(temp, "wb");
  fwrite(components[2], width*height/4, 1, F);
  fclose(F);
  // 这些是后来加上的
  snprintf(temp, 1024, "%s.YUV", filename);
  F = fopen(temp, "wb");
  fwrite(components[0], width, height, F);
  fwrite(components[1], width * height / 4, 1, F);
  fwrite(components[2], width * height / 4, 1, F);
  fclose(F);
}

使用命令行调试：

得到结果如下：

$\text{§}$ 3.2 得到量化矩阵和Huffman表

$\text{§}$ 3.2.1 输出量化矩阵

量化表的代码在build_quantization_table()中，可以添加如下代码后实现输出：

static void build_quantization_table(float *qtable, const unsigned char *ref_table)
{
  /* Taken from libjpeg. Copyright Independent JPEG Group's LLM idct.
   * For float AA&N IDCT method, divisors are equal to quantization
   * coefficients scaled by scalefactor[row]*scalefactor[col], where
   *   scalefactor[0] = 1
   *   scalefactor[k] = cos(k*PI/16) * sqrt(2)    for k=1..7
   * We apply a further scale factor of 8.
   * What's actually stored is 1/divisor so that the inner loop can
   * use a multiplication rather than a division.
   */
  int i, j;
  static const double aanscalefactor[8] = {
     1.0, 1.387039845, 1.306562965, 1.175875602,
     1.0, 0.785694958, 0.541196100, 0.275899379
  };
  const unsigned char *zz = zigzag;

  for (i=0; i<8; i++) {
     for (j=0; j<8; j++) {
       *qtable++ = ref_table[*zz++] * aanscalefactor[i] * aanscalefactor[j];
     }
   }

#if TRACE
  const unsigned char* zz1 = zigzag;
  for (int i = 0; i < 8; i++) {
      for (int j = 0; j < 8; j++) {
          fprintf(p_trace, "%d", ref_table[*zz1++]);
          if (j == 7) {
              fprintf(p_trace, "\n");
          }
      }
  }
#endif
}

得到的结果保存在trace_jpg.txt里面：

$\text{§}$ 3.2.2 输出Huffman码表

构建Huffman的代码在build_huffman_table中，可以直接使用：

$\text{§}$ 3.3 输出某一个DC、AC值图像并统计其概率分布

选择的图片仍然是上述的$1024\times 1024$的图片。

需要在tinyjpg.c中添加代码：

// 保存直流或交流的系数
void save_DC_AC(struct jdec_private* priv, int component)
{
    int DC = (priv->component_infos[component]).DCT[0], AC = priv->component_infos[component].DCT[1];
    DC_buffer[component][DC_pos[component] ++] = DC;
    DC_max[component] = DC_max[component] < DC ? DC : DC_max[component];
    DC_min[component] = DC_min[component] < DC ? DC_min[component] : DC;
    AC_buffer[component][AC_pos[component] ++] = AC;
    AC_max[component] = DC_max[component] < AC ? AC : DC_max[component];
    AC_min[component] = DC_min[component] < AC ? DC_min[component] : AC;
}
// 初始化
void save_Init(struct jdec_private* priv)
{
    for (int i = 0; i < 3; ++i)
    {
        AC_buffer[i] = (int*)malloc(priv->width * priv->height / 16);
        DC_buffer[i] = (int*)malloc(priv->width * priv->height / 16);
        AC_pos[i] = 0; DC_pos[i] = 0;
        AC_max[i] = -1E9; AC_min[i] = 1E9;
        DC_max[i] = -1E9; DC_min[i] = 1E9;

    }
}

同时加入计算概率的函数：

void get_freq()
{
    double* AC_fre[3], * DC_fre[3];
    for (int i = 0; i < 3; ++i)
    {
        AC_fre[i] = (int*)malloc((AC_max[i] - AC_min[i] + 1) * sizeof(double));
        DC_fre[i] = (int*)malloc((DC_max[i] - DC_min[i] + 1) * sizeof(double));
        memset(AC_fre[i], 0, (AC_max[i] - AC_min[i] + 1) * sizeof(double));
        memset(DC_fre[i], 0, (DC_max[i] - DC_min[i] + 1) * sizeof(double));
    }
    for (int i = 0; i < 3; ++i)
    {
        for (int j = 0; j < DC_pos[i]; ++j)
        {
            AC_fre[i][AC_buffer[i][j] - AC_min[i]] ++;
            DC_fre[i][DC_buffer[i][j] - DC_min[i]] ++;
        }

    }
    char AC_file_name[] = "freq_DC_x.csv";
    char DC_file_name[] = "freq_AC_x.csv";
    for (int i = 0; i < 3; ++i)
    {
        DC_file_name[8] = i + '0';
        AC_file_name[8] = i + '0';
        FILE* DC = fopen(DC_file_name, "w+"), * AC = fopen(AC_file_name, "w+");
        fprintf(DC, "value, frequency\n"); fprintf(AC, "value, frequency\n");
        for (int j = 0; j < DC_max[i] - DC_min[i] + 1; ++j)
            fprintf(DC, "%d, %lf\n", j + DC_min[i], (1.0 * DC_fre[i][j]) / DC_pos[i]);
        for (int j = 0; j < AC_max[i] - AC_min[i] + 1; ++j)
            fprintf(AC, "%d, %lf\n", j + AC_min[i], (1.0 * AC_fre[i][j]) / AC_pos[i]);
        fclose(DC); fclose(AC);
    }
}

修改decode_MCU_1x1_3planes

static void decode_MCU_1x1_3planes(struct jdec_private *priv)
{
    // Y
    process_Huffman_data_unit(priv, cY);
#if  TRACE
    save_DC_AC(priv, cY);
#endif //  TRACE
    IDCT(&priv->component_infos[cY], priv->Y, 8);
    // Cb
    process_Huffman_data_unit(priv, cCb);
#if  TRACE
    save_DC_AC(priv, cCb);
#endif //  TRACE
    IDCT(&priv->component_infos[cCb], priv->Cb, 8);

    // Cr
    process_Huffman_data_unit(priv, cCr);
#if  TRACE
    save_DC_AC(priv, cCr);
#endif //  TRACE
    IDCT(&priv->component_infos[cCr], priv->Cr, 8);
}

在tinyjpg.h中加声明：

FILE* DC_coff, DC_freq;
FILE* AC_coff, AC_freq;
FILE* Huffman_code;
FILE* quantization_table;
int* AC_buffer[3], * DC_buffer[3];
int AC_pos[3], DC_pos[3];
int AC_max[3], DC_max[3], AC_min[3], DC_min[3];

得到的结果为：

DC系数图	AC系数图

系数分布：

DC	AC