数据压缩实验：JPEG原理及代码分析

写在前面

首先我想就JPEG说说我自己的认识，然后再结合资料做具体说明。
JPEG编码在我看来其实就是通过舍弃一些人眼识别不太出来的部分，来使得图像没那么精细，然后达到图像的数据量变小的目的，那么我们怎么把这些人眼识别不太出来的部分提取出来呢？
1.我们首先进行色彩变换，这样做的目的是人眼对于亮度信号是很敏感的，而色度信号的稍作改变对于人眼是可以接受的，这就好像人眼在黑暗的条件下颜色的改变看不太出来一样，所以我们要把RBG先变换为YUV，然后我们保留亮度信号Y，对UV分量进行一系列操作。那我们怎么操作呢？
2.我们这里参考了傅里叶的思想，其实傅里叶是将周期信号用不同频率的正余弦基本波来进行分解，那么我们就想了，要是能把图片也进行类似的分解就好了，我们把图像按88分割成一个个小块，挑一个小块出来，我们也可以找到一些基底，然后当然这里我们找了88=64个基底，每个都是一个基本的分量，就好像一个基本的正余弦波一样，我们用这些个基底的系数来表示一张图片，就像下图一样。

这是我从b站视频里面截的，那我们发现，可能我们把这64个基底编上号，也就是说，你这个图片可能是由第1个基地×一个（系数）+第二个基底×一个（系数）+…得来，好，接下来我们怎么才能缩小图片的体积呢？
3.我们发现，这些个颜色变化的特别快的点往往对人眼影响不大，就比如说你把一件白色的羊毛衫放到黑板子上面，这就是说你可能黑白交界处会很明显，但你看不太清羊毛衫上的小孔一样，所以我们可以把那些频率较高的基底赋予一个较低的权重，那这个权重我们怎么统一规定呢？
4.这就要说到量化了，人们针对你想把图片压缩到一个什么水准提出了不同的量化表格，但是他们大多数都有一个规律，就是越靠近频率高的基底他的量化系数就越大，这也就是说，因为我们最终是要除以这个量化表嘛，最终就会和我刚才说的，赋予高频率分量较小的权重达到一个效果。
5.接下来的处理就是，其实我们做完上述步骤，得到的是一个左上角有较多值，然而右下角这一部分大多数哦都被四舍五入赋为零了的这么一个表，那么我们按照Z型给他展开成一个一维数组，就可以按照之前学过的一些基本编码方法进行霍夫曼编码之类的编码，进一步压缩了数据。
6.那么逆过程也是一样的道路反着走就行了，但是由于我们中间有一个这个四舍五入的估算，就是在量化那里，所以我们最终还原出来的图像就和最开始的图像还是有很大差距，视觉感官上看来就是变模糊了。
以上就是我自己关于JPEG编码的全部理解了，以下都是参考这互联网资料写出来的。

JPEG 百度百科简介

JPEG 是Joint Photographic Experts Group(联合图像专家小组)的缩写，是第一个国际图像压缩标准。
JPEG(Joint Photographic Experts Group)是在国际标准化组织(ISO)领导之下制定静态图像压缩标准的委员会，第一套国际静态图像压缩标准ISO 10918-1(JPEG)就是该委员会制定的。由于JPEG优良的品质，使他在短短几年内获得了成功，被广泛应用于互联网和数码相机领域，网站上80%的图像都采用了JPEG压缩标准。

JPEG本身只有描述如何将一个影像转换为字节的数据串流(streaming)，但并没有说明这些字节如何在任何特定的储存媒体上被封存起来。.jpeg/.jpg是最常用的图像文件格式，由一个软件开发联合会组织制定，是一种有损压缩格式，能够将图像压缩在很小的储存空间，图像中重复或不重要的资料会被丢失，因此容易造成图像数据的损伤。尤其是使用过高的压缩比例，将使最终解压缩后恢复的图像质量明显降低，如果追求高品质图像，不宜采用过高压缩比例。但是JPEG压缩技术十分先进，它用有损压缩方式去除冗余的图像数据，在获得极高的压缩率的同时能展现十分丰富生动的图像，换句话说，就是可以用最少的磁盘空间得到较好的图像品质。而且 JPEG是一种很灵活的格式，具有调节图像质量的功能，允许用不同的压缩比例对文件进行压缩，支持多种压缩级别，压缩比率通常在10:1到40:1之间，压缩比越大，品质就越低;相反地，品质就越高。比如可以把1.37Mb的BMP位图文件压缩至20.3KB。当然也可以在图像质量和文件尺寸之间找到平衡点。JPEG格式压缩的主要是高频信息，对色彩的信息保留较好，适合应用于互联网，可减少图像的传输时间，可以支持24bit真彩色，也普遍应用于需要连续色调的图像。

步骤详解

用中，JPEG图像编码算法使用的大多数是离散余弦变换、Huffman编码、顺序编码模式。这样的方式被称之为JPEG基本系统。
1.RGBtoYUV

2.电平偏置下移
这一步因为我们得到的是0-256的电平嘛，我们减去128就让那个中间数为零，就像正弦波进行偏移变换的目的一样，这么做的好处是图像内容平均亮度较高，将0电平移到中间，平均亮度降低， 便于DCT变换量化后直流的系数大大降低，也就降低了数据量。数据形式由无符号数变为有符号数(补码)，单极性数据变为双极性数据。
3.分块
我们进行的DCT是对8*8的小方格进行的，那么我们有一些图片没达到整数，我们就需要去给他填充边缘，填充边缘的方法我们有在电视原理的课程中学习过。
4.DCT
奥妙之处在于，经过DCT，数据中隐藏的规律被发掘了出来,杂乱的数据被转换成几个工整变化的数据。DCT转换后的数组中第一个是一个直线数据，因此又被称为“直流数据”，简称DC，后面的数据被称为“交流数据”，简称AC，这个称呼起源于信号分析中的术语。
在JPEG压缩过程中，经过颜色空间的转换，每一个8X8的图像块，在数据上表现为3个8X8的矩阵，紧接着我们对这三个矩阵做一个二维的DCT转换，二维的DCT转换公式为
正向

反向

这里另外一位博主写的太好了引用一下（引用自jinchao）
离散余弦变换属于傅里叶变换的另外一种形式，没错，就是大名鼎鼎的傅里叶变换。傅里叶是法国著名的数学家和物理学家，1807年，39岁的傅里叶在他的一篇论文里提出了一个想法，他认为任何周期性的函数，都可以分解为为一系列的三角函数的组合，这个想法一开始并没有得到当时科学界的承认，比如当时著名的数学家拉格朗日提出质疑，三角函数无论如何组合，都无法表达带有“尖角”的函数，一直到1822年拉格朗日死后，傅里叶的想法才正式在他的著作《热的解析理论》一书中正式发表。
金子总会闪光，傅里叶变换如今广泛应用于数学、物理、信号处理等等领域，变换除了它在数学上的意义外，还有其哲学上的伟大意义，那就是，世上任何复杂的事物，都可以分解为简单的事物的组合，而这个过程只需要借助数学工具就可以了。但是当年拉格朗日的质疑是正确的，三角函数的确无法表达出尖角形状的函数，不过只要三角函数足够多，可以无限逼近最终结果。比如下面这张动图，就动态描述了一个矩形方波，是如何做傅里叶分析的。

5.量化
我们的问题是，在可以损失一部分精度的情况下，如何用更少的空间存储这些浮点数？
答案是使用量子化（Quantization），简称量化。“量子”这个概念来自于物理学，意思是说连续的能量可以看做是一个个单元体的组合，看起来高端大气，其实很简单，比如游戏中在处理角色面朝方向时，一般并不是使用0到2π这样的浮点数，而是把方向分成16个区间，用0到16这样的整数来表示，这样只用4个bit就足够了。JPEG提供的量子化算法如下。

残差信号经过DCT后，变换系数往往具有较大的动态范围。因此对变换系数进行量化可以有效的减小信号的取值空间，进而获得更好的压缩效果。
量化器主要是利用人眼视觉特性设计而成的矩阵量化DCT系数，减少视觉冗余。
将DCT变换后的临时结果，除以各自量化步长并四舍五入后取整，得到量化系数。
JPEG系统分别规定了亮度分量和色度分量的量化表，色度分量相应的量化步长比亮度分量大。

我们之前从也说过，量化这一步是唯一一步产生误差的一步了。
附上两张标准量化表：


6.将二维矩阵变化为一维的数组
Z字形展开就可以了
在JPEG中，采用曲徊序列，即以矩阵对角线的法线方向作"之"字排列矩阵中的元素。这样做的优点是使得靠近矩阵左上角、值比较大的元素排列在行程的前面，而行程的后面所排列的矩阵元素基本上为0值。
这么做的目的只有一个，就是尽可能把0放在一起，由于0大部分集中在右下角，所以才去这种由左上角到右下角的顺序，经过这种顺序变换，最终矩阵变成一个整数数组。
7.DC系数：差分编码
8*8图像块经过DCT变换后，得到的直流系数有两个特点：数值大、相邻数值变化不大。
根据这个特点，JPEG对DC系数进行差分脉冲调制编码（DPCM）（其中，相邻DC系数之差不在进行量化），对相邻图像块之间的量化DC之差DIFF进行huffman编码。
8.AC系数：游程编码
先将数据之字形扫描，将二位数据扫描成一维数据，定义（run,level）结构，来对此一维数据编码。其中，run表示非0数据前0的个数，level表示非0数据的值。run的最大值为1个，用RRRR表示，level的编码类似DC 系数，先分为16个类别，用4位SSSS表示，再查类内索引（用定长码编码）。即交流huffman码字的前四位表示0的个数，后四位表示该交流分量数值的二进制位数，也即接下来需要读入的位数。

JPEG解码流程

我们首先要分析一下：
JPEG文件结构的分析——段及段的类型
JPEG文件其基本的数据组织方式是段。而段的划分则是由段标记码来划分的，这样的标记码有数十种。且一个以JFIF（JPEG File Interchange Format,JPEG文件交换格式）格式存储的JPEG文件，其段标记码往往以如下的方式出现。

首先是SOI（Start Of Image，0xFFD8），标识一个图像的开始，再便是APP0（application0，0xFFE0）,之后便给出量化表DQT（Define Quantization Table，0xFFDB）然后进入下一层的内容SOF（Start Of Frame，0xFFC0），给出其对应的霍夫曼码表DHT（Define Huffman Table，0xFFC4）最后，进入扫描层SOS（Start Of Scan，0xFFDA），值得注意的是扫描层中为了防止DPCM的差错无限制的传递，还采用了DRI（Define Restart Interval，0xFFDD），来重置差分编码。
而除此之外，JPEG的数据可视为三层，分别为图像层（image），帧层（frame）与扫描（scan）层。
我们知道了这些，就可以大概推断一下解码流程如下：

1.读取文件
2.解析 Segment Marker
解析 SOI
解析 APP0
检查标识“JFIF”及版本并得到一些参数
解析 DQT
得到量化表长度（可能包含多张量化表） 和精度
得到及检查量化表的序号（只能是 0 - 3）
得到量化表内容（64 个数据）
解析 SOF0
得到每个 sample 的比特数、长宽、颜色分量数
得到每个颜色分量的 ID、水平采样因子、垂直采样因子、使用的量化表 序号（与 DQT 中序号对应）
解析 DHT
得到 Huffman 表的类型（AC、DC）、序号 依据数据重建 Huffman 表解析 SOS
得到解析每个颜色分量的 DC、AC 值所使用的 Huffman 表序号（与 DHT 中序号对应）
3.依据每个分量的水平、垂直采样因子计算 MCU 的大小，并得到每个 MCU 中 8*8 宏块的个数
4.对每个 MCU 解码（依照各分量水平、垂直采样因子对 MCU 中每个分量宏块解 码）
①对每个宏块进行 Huffman 解码，得到 DCT 系数
②对每个宏块的 DCT 系数进行 IDCT，得到 Y、Cb、Cr
③遇到 Segment Marker RST 时，清空之前的 DC DCT 系数
5.解析到 EOI，解码结束
6.将 Y、Cb、Cr 转化为需要的色彩空间并保存。

代码分析

main函数
以此实现了如下功能：接受输入输出文件名称进行传递参数，打开TRACEFILE。通过输入参数选择要输出的文件格式，调用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;
}

2.我们接下来就看看他调用的函数
打开输入输出文件，初始化一个结构体，获得文件参数信息，主要调用*tinyjpeg_parse_header函数来解码jpeg图像，最后调用write_yuv这个函数来进行写入yuv文件。
下面文字引用朱轩宇同学的讲解，很详细了，其实我第一遍并没有看懂，有所参考才明白这一段代码的意思。
“首先，将程序加载到内存中然后就关闭；然后解码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;
}

我们可以看到上一段代码中好多指针都指向了一个名字叫做jdec_private的结构体，我们就来看看这个结构体。
3.jdec_private

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

};

在jedc_private 中，定义了指向三个component的指针和三个component结构体，并指示解码过程中的所有信息，量化表，霍夫曼码表以及图像数据。
然后我们再看看这个component.

4.component结构体
其用于参与霍夫曼解码，反量化，IDCT以及彩色空间变换； Hfactor与Vfactor用于说明水平与垂直的采样情况；Q_table指向该component所对应的量化表；AC_table，DC_table用于指向对DC分量与AC分量进行霍夫曼解码时所需的码表；previous_DC是用于保存前一个块的DC值，用于解DPCM；DCT[64]则是对应着一个块的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
};

5.然后这里的霍夫曼表

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

方便直接查表得出结果，省去了大量计算。
6.我们看看他用的函数
（1）.tinyjpeg_parse_header()
解析头文件信息，控制指移动实现核心功能函数parse_JFIF，该函数通过遍历整个文件流来查找0xFF以确定段标识符的位置，并解析该段标识符中的信息，存入jdec结构体中。

/**
* Initialize the tinyjpeg object and prepare the decoding of the stream.
*
* Check if the jpeg can be decoded with this jpeg decoder.
* 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");

//从文件中获取JPEG文件流的信息
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);//分析整个流中的段标识符，并将对应的信息存入priv中

return ret;
}

(2).parse_DQT()
解码量化表：

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

(3).parse_DHT()
解码霍夫曼码表，这一些解码的部分只是编码是的逆过程而已

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

实验结果

其实要输出那个图像，只需要在write_yuv这个函数里面进行修改就可以了，这里我复制代码总有一些问题就直接贴一个图片

输出结果：

输出量化矩阵
量化表的代码在函数build_quantization_table中，添加代码就可以输出了：

实验结果：

霍夫曼码表这个在原来的函数中就能直接输出了：
实验结果：

任务完成