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

实验目的

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

实验原理

JPEG编码原理

鉴于JPEG编码算法可以在提供较大的压缩比的同时，保持较好的显示质量，JPEG逐渐成为最为熟知和广泛使用的数字图像格式和通用标准。                                                                                                                                                                                                                 

 1.零偏置电平下移

对8×8的像块进行零偏置电平下移，即对于灰度级为的像素，减去后将无符号数变为有符号数，减小其取值范围，提高编码效率。

2.8×8 DCT

DCT是一种无损变换，也无法对图像进行压缩

 3.量化

量化是编码流程中唯一会引入误差也是唯一会带来压缩的步骤 

4.DC系数DPCM差分编码

DCT变换后的DC系数特点：系数的值较大；相邻像块的DC系数差值不大（即存在冗余）。根据这个特点，采用DPCM对相邻图像块之间量化DC系数的差值进行编码。 

 5.AC系数Zig-Zag扫描与RLE游程编码

由于DCT后，系数大多数集中在左上角，即低频分量区，因此采用Zig-Zag（之字形）扫描，将系数按频率的高低顺序读出，这样可以出现很多连零的机会，便于进行RLE（Run Length Encoding，游程编码），尤其在最后，如果都是零，给出EOB (End of Block）即可。

 6.Huffman编码

对DC系数DPCM的结果和AC系数RLE的结果进行Huffman编码，共有亮度DC、亮度AC、色差DC、色差AC四张Huffman编码表。

JPEG解码原理

JPEG解码原理完全为编码的逆过程

 JPEG文件格式分析

用Visual Studio二进制编辑器格式打开一个jpeg文件：

1.SOI，Start of Image图像开始：固定标识符0×FFD8

2.EOI，End of Image图像结束：固定标识符0xFFD9 

3.APP0 应用程序保留标记0

4.DQT 定义量化表：固定标识符0xFFDB

重复出现时表示多个量化表（亮度和色度分别量化） 

 5.SOF0一帧图像的开始，固定标识符0xFFC0

长度：0x0011，精度08（每个颜色分量每个像素1byte），图像高度0x0805，图像宽度（0x2506）

6. DHT：哈夫曼表，固定标识符0xFFC4

通常会有4个DHT表（DC系数的亮度、色度；AC系数的亮度、色度） 

001F表长，01高四bit表示表示DC，低四bit为1号表，后续按照huffman标准化编码表示

 实验内容

JPEG解码程序的调试与理解

对JPEG解码程序的理解分为以下三个部分

1.理解三个结构体的设计目的

struct huffman_table：存储Huffman码表。

/* tinyjpeg-internal.h */

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];
};

struct component：储存当前8×8像块中有关解码的信息。

/* tinyjpeg-internal.h */

struct component 
{
  unsigned int Hfactor; // 水平采样因子
  unsigned int Vfactor; // 垂直采样因子
  float* Q_table;   // 指向该8×8块使用的量化表
  struct huffman_table *AC_table;   // 指向该块使用的AC Huffman表
  struct huffman_table *DC_table;   // 指向该块使用的DC Huffman表
  short int previous_DC;    // 前一个块的直流DCT系数
  short int DCT[64];    // DCT系数数组
    
#if SANITY_CHECK
  unsigned int cid;
#endif
};

struct jdec_private：JPEG数据流结构体，用于存储JPEG图像宽高、数据流指针、Huffman码表等内容，并包含struct huffman_table和struct component。

/* tinyjpeg-internal.h */

struct jdec_private
{
  /* Public variables */
  uint8_t *components[COMPONENTS];  /* 分别指向YUV三个分量的三个指针 */
  unsigned int width, height;	/* 图像宽高 */
  unsigned int flags;

  /* Private variables */
  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];		/* 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];
};

2.理解整体框架

从convert_one_image函数切入分析

/* 读取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];

  /* 将JPEG读入缓冲区 */
  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);

tinyjpeg_init函数：创建JPEG数据流结构体，动态分配空间

jdec = tinyjpeg_init();   // 初始化
  if (jdec == NULL)
    exitmessage("Not enough memory to alloc the structure need for decompressing\n");

tinyjpeg_parse_header函数：解析文件头以及其中引用函数各类标识符

 /* 解析JPEG文件头 */
  if (tinyjpeg_parse_header(jdec, buf, length_of_file)<0)
    exitmessage(tinyjpeg_get_errorstring(jdec));

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 marker为起始，否则不是合法的JPEG文件
    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);   // 开始解析JPEG

  return ret;
}

解析各类标识符

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 == 0)
   {
     if (*stream++ != 0xff)
       goto bogus_jpeg_format;
     /* Skip any padding ff byte (this is normal) */
     while (*stream == 0xff)
       stream++;

     marker = *stream++;    // 获取0xFF后的一个字节（即为marker标识符）
     chuck_len = be16_to_cpu(stream);   // length字段
     next_chunck = stream + chuck_len;
     switch (marker)    // 判断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:
	 break;
      }

     stream = next_chunck;  // 解析下一个marker
   }

  if (!dht_marker_found) {
    build_default_huffman_tables(priv);
  }
  return 0;
    
bogus_jpeg_format:
  return -1;
}

重点关注解析DQT

static int parse_DQT(struct jdec_private *priv, const unsigned char *stream)
{
  int qi;   // 量化表ID
  float *table; // 指向量化表
  const unsigned char *dqt_block_end;   // 指向量化表结束位置
  dqt_block_end = stream + be16_to_cpu(stream);
  stream += 2;	// 跳过长度字段

  while (stream < dqt_block_end)	// 检查是否还有量化表
   {
     qi = *stream++;    // 将量化表中系数逐个赋给qi
     table = priv->Q_tables[qi];
     build_quantization_table(table, stream);
     stream += 64;
   }
  return 0;
}

build_quantization_table函数建立量化表

static void build_quantization_table(float *qtable, const unsigned char *ref_table)
{
  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];
     }
   }
}

zigzag数组实现了之字形扫描：

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

解析DHT

获得码表序号以及DC、AC、亮度、色度间的对应关系

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++;
     huff_bits[0] = 0;
     count = 0;
     for (i=1; i<17; i++) {
	huff_bits[i] = *stream++;
	count += huff_bits[i];
     }
#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;
}

获取长、宽信息

/* 计算图像宽高 */
  tinyjpeg_get_size(jdec, &width, &height);

解码

根据上面判断的参数对JPEG数据进行解码，得到解码后的数据，并且对每个宏块进行huffman解码，得到了DCT的系数在进行反DCT

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

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;
  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:   //这种格式使用的decode_mcu_table是decode_mcu_3comp_table
       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;
  }
 
  //mcu的组织，对每个 MCU 解码
  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];   //使用的函数为decode_MCU_1x1_3planes
     convert_to_pixfmt = colorspace_array_conv[0];
  } 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;
 
  //对每个宏块进行 Huffman 解码得到DCT系数
  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;
	    }
	 }
      }
   }
  return 0;
}

写文件

/* 按照指定的输出格式保存输出文件 */
  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;
}

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);
}

得到输出结果：

3.调试TRACE

头文件tinyjpeg.h中：

#define TRACE 1	// 若设为0则可关闭TRACE
#define TRACEFILE "trace_jpeg.txt"	// TRACE文件的文件名

 

 4.输出量化矩阵和Huffman码表

tinyjpeg.h中添加

FILE *DQT_trace;//量化表
FILE *DHT_trace;//huffman表
 
#define DQT_file "DQT.txt"
#define DHT_file "DHT.txt"

tinyjpeg.c中添加

static void build_quantization_table(float *qtable, const unsigned char *ref_table)
{
  ...
  for (i=0; i<8; i++) {
      for (j=0; j<8; j++) {
          /* Added by S.Z.Zheng */
          fprintf(qtabFilePtr, "%-6d", ref_table[*zz]);
          if (j == 7) {
              fprintf(qtabFilePtr, "\n");
          }
          /* Addition ended */
		 ...
      }
  }
  fprintf(qtabFilePtr, "\n\n"); // Added by S.Z.Zheng
}

static int parse_DQT(struct jdec_private *priv, const unsigned char *stream)
{
    ...
    while (stream < dqt_block_end)    // 检查是否还有量化表
    {
        ...
        fprintf(qtabFilePtr, "Quantisation table [%d]:\n", qi);   // 量化表ID
        build_quantization_table(table, stream);
        ...
    }
	...
}

主函数中添加：

int main(int argc, char *argv[])
{
  ...
  /* Added by S.Z.Zheng */
  const char* qtabFileName = "q_table.txt"; // 量化表文件名
  fopen_s(&qtabFilePtr, qtabFileName, "wb");    // 打开文件
  /* Addition ended */
  ...
  fclose(qtabFilePtr);  // Added by S.Z.Zheng
  return 0;
}

输出的量化表如下（0号为DC表，1号为AC表）：

 哈夫曼表如下：

 

 5.输出DC、AC图像及其概率分布

在tinyjpeg_decode函数中修改，DCT[0]即DC图像，DCT[1]~DCT[63]为AC系数：

int tinyjpeg_decode(struct jdec_private *priv, int pixfmt)
{
    FILE* DCFile;///DC系数的存储指针
    FILE* ACFile_1, * ACFile_10, * ACFile_20;AC系数的存储指针
    DCFile = fopen("DC.yuv", "w");
    ACFile_1 = fopen("AC1.yuv", "w");
    ACFile_10 = fopen("AC10.yuv", "w");
    ACFile_20 = fopen("AC10.yuv", "w");
    unsigned char* uvbuf = 128;     //将DC系数和AC系数认为是Y分量，uv分量统一设置为128
    unsigned char* DCbuf, * ACbuf_1, * ACbuf_10, * ACbuf_20;
    int count = 0;    //统计Y分量的数量
    .......
    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);
            //加入DC,AC的数据接入文件
            //DC系数的范围是-512~512，为了使图像能显示，手动+512/4
            DCbuf = (unsigned char)((priv->component_infos->DCT[0] + 512) / 4.0);
 
            fwrite(&DCbuf, 1, 1, DCFile);
            //AC系数，手动+128，调成正值
            ACbuf_1 = (unsigned char)((priv->component_infos->DCT[1] + 128));
            fwrite(&ACbuf_1, 1, 1, ACFile_1);
 
            ACbuf_10 = (unsigned char)((priv->component_infos->DCT[10] + 128));
            fwrite(&ACbuf_10, 1, 1, ACFile_10);
 
            ACbuf_20 = (unsigned char)((priv->component_infos->DCT[20] + 128));
            fwrite(&ACbuf_20, 1, 1, ACFile_20);
	        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;
	        }
	    }
     }
  }
  .......
  //uv分量写进文件
  for (int j = 0; j < count * 0.25 * 2; j++)
  {
      fwrite(&uvbuf, sizeof(unsigned char), 1, DCFile);
      fwrite(&uvbuf, sizeof(unsigned char), 1, ACFile_1);
      fwrite(&uvbuf, sizeof(unsigned char), 1, ACbuf_10);
      fwrite(&uvbuf, sizeof(unsigned char), 1, ACbuf_20);
  }
 
  fclose(DCFile);
  fclose(ACFile_1);
  fclose(ACbuf_10);
  fclose(ACbuf_20);
  return 0;
}

输出结果：

6.计算DC、AC图像的pmf

 

 由图可见，DC系数的波动范围大，分布范围广（具有原图的大部分能量），AC系数能量较小，且集中于128附近。

总结

JPEG编码的问题：

DCT变换以8×8宏块为单位，每个宏块的DCT变换后系数不同、量化不同故产生块效应

当压缩比较高时，压缩相当于经过了一个低通滤波器，时域体现为于sa函数的卷积，卷积的衰减波纹造成像素间的串扰（容易产生边缘模糊等失真现象），在MPEG中体现为“预回声”。