【数据压缩】JPEG原理分析及JPEG解码器的调试

一、实验目的

掌握JPEG编解码系统的基本原理。初步掌握复杂的数据压缩算法实现，并能根据理论分析需要实现所对应数据的输出。

二、实验原理

1、JPEG编解码原理

JPEG编码的过程如上图所示。解码是编码的逆过程。

（1）零偏置电平下移

• 对于灰度级是$2^n$的像素，通过减去$2^{n-1}$，将无符号的整数值变成有符号数。
  对于n=8，即将0-255的值域，通过减去128， 转换为值域在-128~127之间的值。
• 目的：使像素的绝对值出现3位10进制的概率大大减少，提高编码效率。

（2）8×8 DCT变换

对每个单独的彩色图像分量，把整个分量图像分成8×8的图像块，并且这些图象块中的每一个都使用其自己的离散余弦变换单独编码，变换后各分量之间的相关性去除。
正交变换的实质可看作基底函数的分解，如下图所示。

这里是64个余弦函数，通过系数加权，它们可以组合成任何8×8的图像。我们需要做的是计算每个余弦波的系数。左上角全白的是直流系数即DC系数，其它是交流系数即AC系数。从左上角到右下角依次从低频部分到高频部分，由于高频余弦波对图像的贡献不大，通过去除高频部分，使图像空间得到压缩，去除高频部分的过程称为量化。

（3）量化器

• 量化步距是按照系数①所在的位置，②颜色分量 来确定。
• 因为人眼对亮度信号比对色差信号更敏感，因此使用了两种量化表：亮度量化值和色差量化值。
• 根据人眼的视觉特性（对低频敏感，对高频不太敏感）对低频分量采取较细的量化，对高频分量采取较粗的量化。
• 在JPEG标准中，不同的压缩器将根据感觉使用不同的量化表，我们要做的是将每个系数除以相应的量化值，然后四舍五入到最接近的整数。
  基于人的生理感知阈值实验，建议的量化表如下图所示。
  
• 真正的量化表=缩放因子×基本量化表，缩放因子由质量因子决定，质量因子越高越接近原图，压缩率越小。
  质量因子≤ 50：缩放因子= 50 / 质量因子;
  质量因子> 50：缩放因子 = 2 – 质量 因子/ 50

（4）DC系数差分编码

• 8×8图像块经过DCT变换之后得到的DC直流系数有两个特点：
  ➢ 系数的数值比较大
  ➢ 相邻8×8图像块的DC系数值变化不大：冗余

根据这个特点，JPEG算法使用了差分脉冲调制编码(DPCM)技术，对相邻图像块之间量化DC系数的差值DIFF进行编码：$DIF{F}_k=D{C}_k-D{C}_{k-1}$
之后，对DIFF进行Huffman编码。

（5）AC系数游程编码

• 由于经DCT变换后，系数大多数集中在左上角，即低频分量区，因此采用Z字形按频率的高低顺序读出，可以出现很多连零的机会。可以使用游程编码。尤其在最后，如果都是零，给出EOB (End of Block)即可，如下图所示。
  
• 游程编码在JPEG和MPEG编码中规定为：(run, level)
  ➢ 表示连续run个0，后面跟值为level的系数
  ➢ 如：0，2，0，0，3，0，-4，0，0，0，-6，0，0，5，7
  ➢ 表示为(1, 2)， (2, 3) ，…
• 编码：
  ➢ Run: 最多15个，用4位表示RRRR；
  ➢ Level：类似DC，分成16个类别，用4位表示SSSS表示类别号，类内索引；
  ➢ 之后对(RRRR, SSSS)联合用Huffman编码；
  ➢ 对类内索引用定长码编码。

（6）Huffman编码

①DC系数的Huffman编码

对DIFF用Huffman编码：分成类别，类似指数Golomb编码
➢ 类别ID：一元码编码
➢ 类内索引：采用定长码

例：DC=8，上一DC=5，则DIFF=8-5=3，类别ID=2，类内索引=3，则码流=10011。

②AC系数的Huffman编码

对于任何一个RLE的数据对，都可以表示成(RRRR,SSSS)的形式。其中前面的0-16采用自然码RRRR，后面的SSSS则是与DC一致的Huffman分组的编码方式，存储索引。放到码流里的是其组内编码。

• ZRL：表示16个0
• 0，2，0，0，3，0，-4，0，0，0，-6，0，0，5，7,…

2、JPEG文件格式

（1）Segment 的组织形式

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

• 均以 0xFF 开始，后跟 1 byte 的 Marker 和 2 byte 的 Segment length（包含表示Length 本身所占用的 2 byte，不含“0xFF” + “Marker” 所占用的 2 byte）；
• 采用 Motorola 序（相对于 Intel 序），即保存时高位在前，低位在后；
• Data 部分中，0xFF 后若为 0x00，则跳过此字节不予处理；

（2）字段介绍

SOI，Start of Image，图像开始
APP0，Application，应用程序保留标记0
DQT，Define Quantization Table，定义量化表
SOF0，Start of Frame，帧图像开始
DHT，Define Huffman Table，定义哈夫曼表
SOS，Start of Scan，扫描开始 12字节
EOI，End of Image，图像结束 2字节

①SOI & EOI

• SOI ,Start of Image， 图像开始
  标记代码 2字节 固定值0xFFD8
• EOI，End of Image， 图像结束 2字节
  标记代码 2字节 固定值0xFFD9

②APP0 应用程序保留标记0

标记代码 2字节 固定值0xFFE0
包含9个具体字段：

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

③DQT 定义量化表

• 标记代码 2字节 固定值0xFFDB
• 包含9个具体字段：
  ａ．数据长度 2字节 字段ａ和多个字段ｂ的总长度
  ｂ．量化表 数据长度-2字节
  a) 精度及量化表ID 1字节
  高4位：精度，只有两个可选值 0：8位；1：16位　
  低4位：量化表ID，取值范围为0～3
  b) 表项 (64×(精度+1))字节
  例如8位精度的量化表,其表项长度为64×（0+1）=64字节
• 本标记段中，字段ｂ可以重复出现，表示多个量化表，但最多只能出现4次。

④SOF0 帧图像开始

• 标记代码 2字节 固定值0xFFC0
• 包含9个具体字段：
  ａ．数据长度 2字节 ①~⑥六个字段的总长度
  ｂ．精度 1字节 每个数据样本的位数
  通常是8位，一般软件都不支持 12位和16位
  ｃ．图像高度 2字节 图像高度（单位：像素）
  ｄ．图像宽度 2字节 图像宽度（单位：像素）
  ｅ．颜色分量数 1字节 只有3个数值可选
  　　1：灰度图；3：YCrCb或YIQ；4：CMYK
  　　而JFIF中使用YCrCb，故这里颜色分量数恒为3
  ｆ．颜色分量信息 颜色分量数×3字节（通常为9字节）
  a)颜色分量ID 1字节
  b)水平/垂直采样因子 1字节
  高4位：水平采样因子　低4位：垂直采样因子
  c) 量化表 1字节 当前分量使用的量化表的ID

⑤DHT，定义哈夫曼表

• 标记代码 2字节 固定值0xFFC4

• 包含2个具体字段：
  ａ． 数据长度 2字节
  ｂ．huffman表 数据长度-2字节

• 表ID和表类型 1字节
  高4位：类型，只有两个值可选　0：DC直流；1：AC交流
  低4位：哈夫曼表ID，注意，DC表和AC表分开编码

• 不同位数的码字数量 16字节

• 编码内容 16个不同位数的码字数量之和（字节）
  本标记段中，字段ｂ可以重复出现（一般4次），也可以只出现1次。

⑥SOS，扫描开始

• 标记代码 2字节 固定值0xFFDA

• 包含2个具体字段：
  ａ．数据长度 2字节 ａ~ｄ两个字段的总长度
  ｂ．颜色分量数 1字节 应该和SOF中的字段⑤的值相同，即：1：灰度图是；3： YCrCb或YIQ；4：CMYK。
  ｃ．颜色分量信息
  a) 颜色分量ID 1字节
  b) 直流/交流系数表号 1字节
  高4位：直流分量使用的哈夫曼树编号
  低4位：交流分量使用的哈夫曼树编号
  ｄ．压缩图像数据
  a)谱选择开始 1字节 固定值0x00
  b)谱选择结束 1字节 固定值0x3F
  c)谱选择 1字节 在基本JPEG中总为00

3、JPEG解码流程

（1）读取文件
（2）解析segmentmarker：依次解析出SOI、APP0、DCT、SCF0、DHT、SOS、EOI
解析SOI
解析APP0

• 检查标识“JFIF”及版本
• 得到一些参数

解析DCT

• 得到量化表长度（可能包含多张量化表）
• 得到量化表的精度
• 得到及检查量化表的序号（只能是 0-3）
• 得到量化表内容（64 个数据）

解析SCF0

• 得到每个 sample 的比特数、长宽、颜色分量数
• 得到每个颜色分量的 ID、水平采样因子、垂直采样因子、使用的量化表序号（与 DQT 中序号对应）
• 解析 DHT

解析DHT

• 得到 Huffman 表的类型（AC、DC）、序号
• 依据数据重建 Huffman 表

解析SOS

• 得到解析每个颜色分量的 DC、AC 值所使用的 Huffman 表序号（与 DHT 中序号对应）

（3）依据每个分量的水平采样、垂直采样因子计算MCU（微控制单元，即RAM，即要开辟的图像内存）的大小，并得到每个MCU中８＊８宏块的个数。
（４）对每个MCU解码：对每个宏块进行huffman解码得到DCT系数，进行IDCT得到变换之前的数据。
（５）解析到EOI，解码结束
（６）将得到的Ｙ、Cb、Cr转换成需要的色彩空间并保存。

三、实验过程

逐步调试JPEG解码器程序。将输入的JPG文件进行解码，将输出文件保存为可供YUVViewer观看的YUV文件。
输入文件如下。

工程中有三个头文件：stdint.h、tinypeg.h、tinyjpeg-internal.h以及三个.c文件：jidctflc.c、loadjpeg.c、tinyjpeg.c。

1、关键代码分析

（1）tinyjpeg-internal.h中的结构体

①struct huffman_table

创建一个快速查找表用于快速解码，如果查找失败则用慢速查找表。目的是加快解码过程。

struct huffman_table
{
  /*快速查表，用HUFFMAN_HASH_SIZE个比特可以快速找到符号，如果符号<0，那么需要查看慢速查找表 */
  short int lookup[HUFFMAN_HASH_SIZE];
  /* 给出码字长度，输出权值对应的码长 */
  unsigned char code_size[HUFFMAN_HASH_SIZE];
  /*给出在慢速查找表中没有出现的码字，计算256值是否足以存储所有值*/
  uint16_t slowtable[16-HUFFMAN_HASH_NBITS][256];
};

②struct component

用于存储一个MCU块的信息，每处理一个新的MCU块后，信息都会更新。用它来保存一些解码过程中需要用到的信息，如量化表、前一个直流的值和当前MCU的DCT块从而使解码更加简单。

struct component 
{
  unsigned int Hfactor; /* 水平采样情况 */
  unsigned int Vfactor; /* 垂直采样情况 */
  float *Q_table;		/* 指向量化表 */
  struct huffman_table *AC_table; /* 指向直流huffman_table */
  struct huffman_table *DC_table; /* 指向交流huffman_table */
  short int previous_DC;	/* 前一个直流系数 */
  short int DCT[64];		/* DCT系数 */
#if SANITY_CHECK
  unsigned int cid;
#endif
};

③struct jdec_private

该结构体是解码过程中频繁用到的东西，起到整合基本信息、连接各个部分的作用。它用于存储图像的基本信息和各通道信息。包含了前两个结构体的内容。由于是采用指针方式，并不占用额外的空间。

struct jdec_private
{
  /* 公有变量 */
  uint8_t *components[COMPONENTS];
  unsigned int width, height;	/* 图像大小 */
  unsigned int flags;

  /* 私有变量 */
  const unsigned char *stream_begin, *stream_end;
  unsigned int stream_length;

  const unsigned char *stream;	/* 指向当前流的指针 */
  unsigned int reservoir, nbits_in_reservoir;

  struct component component_infos[COMPONENTS];	/* 各通道信息 */
  float Q_tables[COMPONENTS][64];		/* 三通道量化表 */
  struct huffman_table HTDC[HUFFMAN_TABLES];	/* 直流huffman tables   */
  struct huffman_table HTAC[HUFFMAN_TABLES];	/* 交流huffman tables   */
  int default_huffman_table_initialized;
  int restart_interval;
  int restarts_to_go;				/* 剩余MCU */
  int last_rst_marker_seen;			/* Rst marker每次递增 */

  /* IDCT之后用于存储每个component的临时空间  */
  uint8_t Y[64*4], Cr[64], Cb[64];

  jmp_buf jump_state;
  /* 内部指针用于颜色空间转换，请勿修改!!! */
  uint8_t *plane[COMPONENTS];

};

（2）解析函数

①convert_one_image( )

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

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://输出格式为yuv420
      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;
}

②tinyjpeg_parse_header( )

控制指针移动，解析JPEG头文件，调用parse_JFIF()

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))//jpeg文件以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;
}

③parse_JFIF( )

在到达SOS之前，循环调用parse_SOF()、parse_DQT()、parse_SOS()、parse_DHT()、parse_DRI()。

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 */
  while (!sos_marker_found)
   {
     if (*stream++ != 0xff)//获取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;
}

④parse_*( )

parse_SOF()、parse_DQT()、parse_SOS()、parse_DHT()、parse_DRI()分别用于解析各字段。

（3）计算MCU

在此函数中，对每个块的水平和垂直采样情况进行解析，计算MCU。同时对不同MCU进行处理。

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) {//根据不同的输出格式确定MCU
     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;//初始化MCU宽高8px（4：4：4）
  //若水平垂直采样因子均为1
  if ((priv->component_infos[cY].Hfactor | priv->component_infos[cY].Vfactor) == 1) {
     decode_MCU = decode_mcu_table[0];//MCU包含1个Y
     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) {//若水平采样因子为1、垂直采样因子为2
     decode_MCU = decode_mcu_table[1];//MCU包含2个Y
     convert_to_pixfmt = colorspace_array_conv[1];
     ystride_by_mcu = 16;//MCU高16px、宽8px
#if TRACE
     fprintf(p_trace,"Use decode 1x2 sampling (not supported)\n");
	 fflush(p_trace);
#endif
  } else if (priv->component_infos[cY].Vfactor == 2) {//若水平垂直采样因子均为2
     decode_MCU = decode_mcu_table[3];//MCU包含4个Y
     convert_to_pixfmt = colorspace_array_conv[3];
     xstride_by_mcu = 16;//宽高均为16px
     ystride_by_mcu = 16;
#if TRACE 
	 fprintf(p_trace,"Use decode 2x2 sampling\n");
	 fflush(p_trace);
#endif
  } else {//若水平采样因子为2、垂直采样因子为1
     decode_MCU = decode_mcu_table[2];//MCU包含2个Y
     convert_to_pixfmt = colorspace_array_conv[2];
     xstride_by_mcu = 16;//MCU宽16px、高8px
#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;

  /* 对每个像块进行解码 (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;
}

（4）MCU解码

对每个MCU进行解码（依照各分量水平、垂直采样因子对 MCU中每个分量宏块解码）
①对每个宏块进行 Huffman解码，得到DCT系数
②对每个宏块的 DCT系数进行 IDCT，得到 Y、Cb 、Cr
③遇到 Segment Marker RST时，清空之前的DC 、DCT系数
以decode_MCU_2x1_1plane( )为例。

①decode_MCU_2x1_1plane( )

/*
 * Decode a 2x1
 *  .-------.
 *  | 1 | 2 |
 *  `-------'
 */
static void decode_MCU_2x1_1plane(struct jdec_private *priv)
{
  // Y
  process_Huffman_data_unit(priv, cY);
  IDCT(&priv->component_infos[cY], priv->Y, 16);//块1
  process_Huffman_data_unit(priv, cY);
  IDCT(&priv->component_infos[cY], priv->Y+8, 16);//块2
  // Cb
  process_Huffman_data_unit(priv, cCb);
  // Cr
  process_Huffman_data_unit(priv, cCr);
}

②process_Huffman_data_unit( )

对包含DCT系数的单个块进行Huffman解码。

/**
 *
 * Decode a single block that contains the DCT coefficients.
 * The table coefficients is already dezigzaged at the end of the operation.
 *
 */
static void process_Huffman_data_unit(struct jdec_private *priv, int component)
{
  unsigned char j;
  unsigned int huff_code;
  unsigned char size_val, count_0;

  struct component *c = &priv->component_infos[component];
  short int DCT[64];


  /* Initialize the DCT coef table */
  memset(DCT, 0, sizeof(DCT));

  /* DC coefficient decoding */
  huff_code = get_next_huffman_code(priv, c->DC_table);
  //trace("+ %x\n", huff_code);
  if (huff_code) {
  	 //解码得到当前块与前一块DC系数的差值
     get_nbits(priv->reservoir, priv->nbits_in_reservoir, priv->stream, huff_code, DCT[0]);
     DCT[0] += c->previous_DC;//差值与前一DC系数相加，得到当前块DC系数
     c->previous_DC = DCT[0];//更新前一DC系数值
  } else {
     DCT[0] = c->previous_DC;
  }

  /* AC coefficient decoding */
  j = 1;
  while (j<64)
   {
     huff_code = get_next_huffman_code(priv, c->AC_table);
     //trace("- %x\n", huff_code);

     size_val = huff_code & 0xF;
     count_0 = huff_code >> 4;

     if (size_val == 0)
      { /* RLE */
	if (count_0 == 0)
	  break;	/* EOB found, go out */
	else if (count_0 == 0xF)
	  j += 16;	/* skip 16 zeros */
      }
     else
      {
	j += count_0;	/* skip count_0 zeroes */
	if (__unlikely(j >= 64))
	 {
	   snprintf(error_string, sizeof(error_string), "Bad huffman data (buffer overflow)");
	   break;
	 }
	get_nbits(priv->reservoir, priv->nbits_in_reservoir, priv->stream, size_val, DCT[j]);
	j++;
      }
   }

  for (j = 0; j < 64; j++)
    c->DCT[j] = DCT[zigzag[j]];
}

（5）色彩空间转换

以YCrCb -> YUV420P (2x1)为例。

/**
 *  YCrCb -> YUV420P (2x1)
 *  .-------.
 *  | 1 | 2 |
 *  `-------'
 */
static void YCrCB_to_YUV420P_2x1(struct jdec_private *priv)
{
  unsigned char *p;
  const unsigned char *s, *y1;
  unsigned int i;

  p = priv->plane[0];
  y1 = priv->Y;
  for (i=0; i<8; i++)
   {
     memcpy(p, y1, 16);
     p += priv->width;
     y1 += 16;
   }

  p = priv->plane[1];
  s = priv->Cb;
  for (i=0; i<8; i+=2)
   {
     memcpy(p, s, 8);
     s += 16; /* Skip one line */
     p += priv->width/2;
   }

  p = priv->plane[2];
  s = priv->Cr;
  for (i=0; i<8; i+=2)
   {
     memcpy(p, s, 8);
     s += 16; /* Skip one line */
     p += priv->width/2;
   }
}

2、生成YUV图像

在保存yuv文件的函数write_yuv( )中末尾处添加如下代码：

  snprintf(temp, 1024, "%s.YUV", filename);
  F = fopen(temp, "wb");
  fwrite(components[0], width, height, F);//写入Y
  fwrite(components[1], width * height / 4, 1, F);//写入U
  fwrite(components[2], width * height / 4, 1, F);//写入V
  fclose(F);

生成如下文件：

用YUVViewer打开output.YUV，如下图所示：

3、TRACE的目的和含义

trace定义在tinyjpeg.h末尾处。
trace的目的主要是输出中间变量或者错误信息，有利于方便高效地理解程序运行的过程及排查错误。

（1）trace的打开和关闭

#define TRACE 1//0为关闭，1为开启
#define  TRACEFILE "trace_jpeg.txt"//输出文件名

（2）各TRACE位置及输出内容

若需输出中间变量或错误信息，仅需在对应位置按照需求修改对应内容即可。

#if TRACE
  /**/
  fflush(p_trace);
#endif

①以txt文件输出所有的量化矩阵

在tinyjpeg.h中添加如下声明：

FILE* quantization_table;

在主程序中为量化矩阵分配文件指针。

#if TRACE
  quantization_table = fopen("quantization_table.txt", "w+");
  if (quantization_table == NULL)
  {
      printf("quantization_table file open error!");
  }
#endif

在parse_DQT( )添加如下代码：

static int parse_DQT(struct jdec_private* priv, const unsigned char* 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
#if TRACE
		fprintf(quantization_table, "original quantization table [%d]\n", qi);
#endif
		table = priv->Q_tables[qi];
		build_quantization_table(table, stream);
		stream += 64;
#if TRACE
		fprintf(quantization_table, "scaled quantization table [%d]\n", qi);
		for (int i = 0; i < 8; ++i)
			for (int j = 0; j < 8; ++j)
				if (j == 7) fprintf(quantization_table, "%10.3f\n", table[i * 8 + j]);
				else fprintf(quantization_table, "%10.3f\t", table[i * 8 + j]);
#endif
	...
	}

在build_quantization_table( )中添加如下代码：

static void build_quantization_table(float* qtable, const unsigned char* ref_table)
{
	...
	for (i = 0; i < 8; i++) {
		for (j = 0; j < 8; j++) {
			*qtable++ = ref_table[*zz++] * aanscalefactor[i] * aanscalefactor[j];
		}
	}
#if TRACE
	for (int i = 0; i < 8; ++i)
		for (int j = 0; j < 8; ++j)
			if (j == 7) fprintf(quantization_table, "%10d\n", (int)ref_table[zigzag[i * 8 + j]]);
			else fprintf(quantization_table, "%10d\t", (int)ref_table[zigzag[i * 8 + j]]);
#endif
}

结果输出为文件quantization_table.txt：

②以txt文件输出所有的HUFFMAN码表

在parse_DHT中已有如下代码：

static int parse_DHT(struct jdec_private* priv, const unsigned char* stream)
{
	...
#if TRACE
		fprintf(p_trace, "Huffman table %s[%d] length=%d\n", (index & 0xf0) ? "AC" : "DC", index & 0xf, count);
		fflush(p_trace);
#endif
	...
}

在build_huffman_table中已有如下代码：

static void build_huffman_table(const unsigned char* bits, const unsigned char* vals, struct huffman_table* table)
{
	...
#if TRACE
		fprintf(p_trace, "val=%2.2x \tcode=%8.8x \tcodesize=%2.2d\n", val, code, code_size);
		fflush(p_trace);
#endif
	...
}

结果输出在trace_jpeg.txt，同时依照quantization_table.txt方式将结果存放在新建文件中，只需将p_trace修改为新文件指针。由于截图过长，以下将内容展示。

• 直流码表0：

Huffman table DC[0] length=10
val=04 	code=00000000 	codesize=02
val=05 	code=00000001 	codesize=02
val=06 	code=00000002 	codesize=02
val=03 	code=00000006 	codesize=03
val=02 	code=0000000e 	codesize=04
val=01 	code=0000001e 	codesize=05
val=00 	code=0000003e 	codesize=06
val=09 	code=0000007e 	codesize=07
val=07 	code=000000fe 	codesize=08
val=08 	code=000001fe 	codesize=09

• 直流码表1：

Huffman table DC[1] length=11
val=04 	code=00000000 	codesize=02
val=05 	code=00000001 	codesize=02
val=06 	code=00000002 	codesize=02
val=03 	code=00000006 	codesize=03
val=02 	code=0000000e 	codesize=04
val=01 	code=0000001e 	codesize=05
val=00 	code=0000003e 	codesize=06
val=07 	code=0000007e 	codesize=07
val=0a 	code=000000fe 	codesize=08
val=09 	code=000001fe 	codesize=09
val=08 	code=000003fe 	codesize=10

• 交流码表0：

Huffman table AC[0] length=43
val=00 	code=00000000 	codesize=02
val=01 	code=00000002 	codesize=03
val=03 	code=00000003 	codesize=03
val=02 	code=00000008 	codesize=04
val=04 	code=00000009 	codesize=04
val=05 	code=0000000a 	codesize=04
val=11 	code=0000000b 	codesize=04
val=21 	code=0000000c 	codesize=04
val=22 	code=0000001a 	codesize=05
val=31 	code=0000001b 	codesize=05
val=61 	code=0000001c 	codesize=05
val=06 	code=0000003a 	codesize=06
val=12 	code=0000003b 	codesize=06
val=a1 	code=0000003c 	codesize=06
val=32 	code=0000007a 	codesize=07
val=41 	code=0000007b 	codesize=07
val=62 	code=0000007c 	codesize=07
val=13 	code=000000fa 	codesize=08
val=51 	code=000000fb 	codesize=08
val=23 	code=000001f8 	codesize=09
val=42 	code=000001f9 	codesize=09
val=71 	code=000001fa 	codesize=09
val=81 	code=000001fb 	codesize=09
val=91 	code=000001fc 	codesize=09
val=15 	code=000003fa 	codesize=10
val=52 	code=000003fb 	codesize=10
val=63 	code=000003fc 	codesize=10
val=07 	code=000007fa 	codesize=11
val=14 	code=000007fb 	codesize=11
val=33 	code=000007fc 	codesize=11
val=53 	code=000007fd 	codesize=11
val=16 	code=00000ffc 	codesize=12
val=43 	code=00000ffd 	codesize=12
val=08 	code=00001ffc 	codesize=13
val=b1 	code=00001ffd 	codesize=13
val=34 	code=00003ffc 	codesize=14
val=c1 	code=00003ffd 	codesize=14
val=24 	code=00007ffc 	codesize=15
val=d1 	code=0000fffa 	codesize=16
val=09 	code=0000fffb 	codesize=16
val=72 	code=0000fffc 	codesize=16
val=f0 	code=0000fffd 	codesize=16
val=a2 	code=0000fffe 	codesize=16

• 交流码表1

Huffman table AC[1] length=28
val=00 	code=00000000 	codesize=02
val=04 	code=00000001 	codesize=02
val=01 	code=00000004 	codesize=03
val=02 	code=00000005 	codesize=03
val=03 	code=00000006 	codesize=03
val=31 	code=0000001c 	codesize=05
val=61 	code=0000001d 	codesize=05
val=11 	code=0000003c 	codesize=06
val=12 	code=0000003d 	codesize=06
val=05 	code=0000007c 	codesize=07
val=21 	code=0000007d 	codesize=07
val=13 	code=000000fc 	codesize=08
val=14 	code=000000fd 	codesize=08
val=41 	code=000000fe 	codesize=08
val=51 	code=000001fe 	codesize=09
val=06 	code=000007fc 	codesize=11
val=22 	code=000007fd 	codesize=11
val=32 	code=000007fe 	codesize=11
val=07 	code=00001ffc 	codesize=13
val=15 	code=00001ffd 	codesize=13
val=42 	code=00001ffe 	codesize=13
val=08 	code=0000fff8 	codesize=16
val=71 	code=0000fff9 	codesize=16
val=23 	code=0000fffa 	codesize=16
val=24 	code=0000fffb 	codesize=16
val=33 	code=0000fffc 	codesize=16
val=81 	code=0000fffd 	codesize=16
val=a1 	code=0000fffe 	codesize=16

4、输出DC图像和某个AC图像并统计其概率分布

（1）输出DC、AC图像

在得到的DCT数组中存储了图像DCT变换后的数据，将这些数据输出即可，其中DCT[0]为图像的直流分量，其余为交流分量，实验中取DCT[1]作为输出的AC分量值。
对于8x8的宏块，取每一个宏块DCT变换后左上角的第一个值为DC分量，输出图像的大小为128x128（原图像为1024x1024）。
为了便于观察，将DC分量的值（范围-512-512）调整到 0~255。AC分量的值较小，+128处理。
在程序中添加输出DC、AC图像的代码。
在tinyjpeg.h中添加全局变量：

FILE *DCfile;
FILE *ACfile;

在main函数中添加如下代码：

DCfile = fopen("dc.yuv","w");
ACfile = fopen("ac.yuv", "w");
........
fclose(DCfile);
fclose(ACfile);

在tinyjpeg_decode函数中添加如下代码：

int tinyjpeg_decode(struct jdec_private *priv, int pixfmt)
{
	unsigned char *DCbuf, *ACbuf; 
	unsigned char *uvbuf = 128;
	int count=0;
	...
	decode_MCU(priv);

	DCbuf= (unsigned char)((priv->component_infos->DCT[0] + 512.0) / 4 + 0.5);  
	ACbuf= (unsigned char)(priv->component_infos->DCT[1] + 128);   
	fwrite(&DCbuf, 1, 1,DCfile);
	fwrite(&ACbuf, 1, 1,ACfile);
	count++;
	....
 	for (int j = 0; j < count *0.25 * 2; j++)
  	{
	  	fwrite(&uvbuf, sizeof(unsigned char), 1, DCfile);
	  	fwrite(&uvbuf, sizeof(unsigned char), 1, ACfile);
  	}

 	return 0;
}

用YUVViewer打开dc.yuv和ac.yuv，如下图所示：

（2）DC、AC图像的概率分布图

DC和AC的概率分布图如下所示：

python代码如下所示：

import numpy as np                
import matplotlib.pyplot as plt   
import math                      

f = open(r"dc.yuv", "rb")
data = f.read()
f.close()
data = [int(x) for x in data]
data_Y=[]

for i in range(0,128*128):
    data_Y.append(data[i])

count_Y=np.zeros(255)
for i in data_Y:
    count_Y[i-1]=count_Y[i-1]+1

for i in range(0,255):
    count_Y[i] = count_Y[i] / (128*128)

plt.plot(count_Y)
plt.title("DC Frequency")
plt.grid()
plt.show()

四、实验问题

1、问题：PDB格式不兼容

解决方法：点击 生成-清理解决方案。

---

参考

JPEG编解码原理及C++调试
JPEG原理分析及JPEG解码器的解析
JPEG编码原理与解码分析
利用python分析RGB图像及YUV图像三个通道的概率分布并求其熵