【实验4】 DPCM 压缩系统的实现和分析

---

一、实验目的

掌握DPCM编解码系统的基本原理。初步掌握实验用C/C++/Python等语言编程实现DPCM编码器，并分析其压缩效率。

二、实验原理和实验要求

1.DPCM编解码原理

  DPCM是差分预测编码调制的缩写，是比较典型的预测编码系统。在DPCM系统中，需要注意的是预测器的输入是已经解码以后的样本。之所以不用原始样本来做预测，是因为在解码端无法得到原始样本，只能得到存在误差的样本。因此，在DPCM编码器中实际内嵌了一个解码器，如编码器中虚线框中所示。
  在一个DPCM系统中，有两个因素需要设计：预测器和量化器。理想情况下，预测器和量化器应进行联合优化。实际中，采用一种次优的设计方法：分别进行线性预测器和量化器的优化设计。

2.DPCM编码系统的设计

• 在本次实验中，采用固定预测器和均匀量化器。
• 本实验的目标是验证DPCM编码的编码效率。
    首先读取一个256级的灰度图像，采用自己设定的预测方法计算预测差，并对预测误差进行8比特均匀量化(基本要求)。还可对预测误差进行1比特、2比特和4比特的量化设计(提高要求)。在DPCM编码器实现的过程中可同时输出预测误差图像和重建图像。
    将预测误差图像写入文件并将该文件输入Huffman编码器，得到输出码流、给出概率分布图并计算压缩比。将原始图像文件输入Huffman编码器，得到输出码流、给出概率分布图并计算压缩比。
    最后比较两种系统（1.DPCM+熵编码和2.仅进行熵编码）之间的编码效率（压缩比和图像质量）。压缩质量以PSNR进行计算。

三、实验关键代码

1.预测+量化

• 使用左侧预测的方法

void DPCM(unsigned char* ori_yBuf, unsigned char* pre_yBuf, unsigned char* rec_yBuf, int width, int height, int bitdepth)
{
	double num = pow(2, (int)(9 - bitdepth));
	for (int i = 0; i < height; i++)    //每一行
	{
		for (int j = 0; j < width; j++)  //第一列
		{
			if (j == 0) //使用左侧预测的方法，第一列的像素没有参考值
			{
				//假设预测值为128
				//计算预测误差，对预测误差进行量化
				pre_yBuf[i * width] = (unsigned char)(((ori_yBuf[i * width] - 128) + 255) / num);
				//计算重建值，重建值=反量化后的误差+预测值
				rec_yBuf[i * width] = (unsigned char)((pre_yBuf[i * width] - 255 / num) * num + 128);

				
			}
			else
			{   //选取前一像素的重建值作为预测值
				//计算预测误差，对预测误差进行量化
				pre_yBuf[i * width + j] = (unsigned char)(((ori_yBuf[i * width + j] - rec_yBuf[i * width + j - 1]) + 255) / num);
				//计算重建值，重建值=反量化后的误差+预测值
				rec_yBuf[i * width + j] = (unsigned char)((pre_yBuf[i * width + j]- 255/num) * num + rec_yBuf[i * width + j - 1]);
				
			}
			//防止溢出
			pre_yBuf[i * width + j] = (unsigned char)(pre_yBuf[i * width + j] * num / 2);
			if (rec_yBuf[i * width + j] > 255)
				rec_yBuf[i * width + j] = 255;
			if (rec_yBuf[i * width + j] < 0)
				rec_yBuf[i * width + j] = 0;
		
		}
	}
	
}

2.PSNR计算压缩质量

PSNR(Peak Signal to Noise Ratio),峰值信噪比

• PSNR是一种评价图像的的客观标准，其公式为：

$$PSNR=10\ast {\lg }_{}(\frac{MA{X}^2}{MSE})$$
    其中：MAX指最大值，例如在8bit灰度图中MAX=255;
       MSE指均方误差，其公式为：
$$MSE=\frac{{\sum }_{i=0}^M{\sum }_{j=0}^N[f(i,j)-f^{\prime }(i,j)]}{MN}$$
      在本实验中，f(x,y)和f’(x,y)分别指原图像和重建图像，M和N分别指图像的高和宽。

• 一般来说，PSNR值在20-40dB左右，值越大图像还原质量越好。

void PSNR(unsigned char* ori_yBuf, unsigned char* rec_yBuf, int width ,int height){
	double psnr = 0, MSE = 0;
	for (int i = 0; i < height; i++)
	{
		for (int j = 0;j < width; j++) {
			MSE += (ori_yBuf[i * width + j] - rec_yBuf[i * width + j]) * (ori_yBuf[i * width + j] - rec_yBuf[i * width + j]);
		}
	}
	MSE = (double)MSE/ (double)(height * width);
	int MAX = 255; //本次实验原始图像都是8bit
	psnr = 10 * log10((double)(MAX*MAX)/ MSE);
	cout << "psnr=" << psnr << endl;
}

四、实验步骤与运行结果

• 实验图像准备

    实验使用6张bmp图片（256*256,8bit）
    使用之前实验编写的bmp2yuv代码将图像转换成yuv(4:2:0)格式

1、 DPCM

• 在命令行参数中进行设置
  从左到右依次为原始图像名称，预测误差图像名称，重建图像名称，图像宽度，图像高度，量化bit数。
  
• 运行结果（预测误差选择8bit均匀量化）

• 对比不同bit量化

2、比较图像质量（PSNR）

图像	8bit	4bit	2bit	1bit
Camman256B	27.0756	9.91553	7.19489	7.05763
Clown256B	51.1614	19.4272	9.3693	8.96962
Fruit256B	44.1822	17.7163	8.77114	7.62002
Noise256B	14.8862	-	-	-
Odie256B	18.2415	-	-	-
Zone256B	17.0843	-	-	-

3、Huffman熵编码

    将原始图像文件和预测误差图像输入Huffman编码器，得到输出码流。

- 概率分布图

• 图像Camman
  
  
• 图像Clown
  
  
• 图像Fruit
  
  
• 图像Odie
  
  
• 图像Noise
  
  
• 图像Zone
  
  

- 压缩比

由上表可知，图像Camman、Clown、Fruit、Odie的DPCM+huffman编码的压缩比，比只进行huffman编码的压缩比要大得多。而图像Noise和Zone则相反，DPCM+huffman的压缩比要小于等于只进行huffman编码的压缩比。

分析图像符号的概率分布图可知，图像Camman、Clown、Fruit、Odie概率分布图中，像素值出现概率大的相对集中，即像素间相关性比较强，因此这些图经过DPCM后的概率分布后的概率分布图，符号都集中在127附近，符号种类数大大减少，出现的概率大大增加。而huffman编码系统就喜欢这种分布不均匀的信源，因此编码效率较高。
而Noise、Zone的信源概率分布比较分散，经huffman编码后，压缩效果并不明显，由于符号间相关性不大，经DPCM之后，信源符号个数并没有减少，反而增加了（概率分布图变密），因此经过DPCM+huffman编码后，压缩效果不是很好，比单独进行huffman编码还要低。

---

代码

#include <iostream>
#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>
#include<math.h>
#include<stdint.h>
using namespace std;


//预测+量化
void DPCM(unsigned char* ori_yBuf, unsigned char* pre_yBuf, unsigned char* rec_yBuf, int width, int height, int bitdepth)
{
	double num = pow(2, (int)(9 - bitdepth));
	for (int i = 0; i < height; i++)    //每一行
	{
		for (int j = 0; j < width; j++)  //第一列
		{
			if (j == 0) //使用左侧预测的方法，第一列的像素没有参考值
			{
				//假设预测值为128
				//计算预测误差，对预测误差进行量化
				pre_yBuf[i * width] = (unsigned char)(((ori_yBuf[i * width] - 128) + 255) / num);
				//计算重建值，重建值=反量化后的误差+预测值
				rec_yBuf[i * width] = (unsigned char)((pre_yBuf[i * width] - 255 / num) * num + 128);

				
			}
			else
			{   //选取前一像素的重建值作为预测值
				//计算预测误差，对预测误差进行量化
				pre_yBuf[i * width + j] = (unsigned char)(((ori_yBuf[i * width + j] - rec_yBuf[i * width + j - 1]) + 255) / num);
				//计算重建值，重建值=反量化后的误差+预测值
				rec_yBuf[i * width + j] = (unsigned char)((pre_yBuf[i * width + j]- 255/num) * num + rec_yBuf[i * width + j - 1]);
				
			}
			//防止溢出
			pre_yBuf[i * width + j] = (unsigned char)(pre_yBuf[i * width + j] * num / 2);
			if (rec_yBuf[i * width + j] > 255)
				rec_yBuf[i * width + j] = 255;
			if (rec_yBuf[i * width + j] < 0)
				rec_yBuf[i * width + j] = 0;
		
		}
	}
	
}
//计算压缩质量
void PSNR(unsigned char* ori_yBuf, unsigned char* rec_yBuf, int width ,int height){
	double psnr = 0, MSE = 0;
	for (int i = 0; i < height; i++)
	{
		for (int j = 0;j < width; j++) {
			MSE += (ori_yBuf[i * width + j] - rec_yBuf[i * width + j]) * (ori_yBuf[i * width + j] - rec_yBuf[i * width + j]);
		}
	}
	MSE = MSE/(height * width);
	double MAX = 255;
	psnr = 10 * log10((MAX*MAX)/ MSE);
	cout << "psnr=" << psnr << endl;
}

//得到图像的灰度值分布
void GetFrequency(unsigned char* buffer, double* frequency, int height, int width)
{
	int length = height * width;
	for (int i = 0; i < length; i++) frequency[buffer[i]] ++;
	for (int i = 0; i < 256; i++) frequency[i] /= length;
}


int main(int argc, char* argv[])
{
	char* ori_yuvfilename = NULL;
	char* pre_yuvfilename = NULL;
	char* rec_yuvfilename = NULL;

	FILE* ori_y_file = NULL;
	FILE* pre_y_file = NULL;
	FILE* rec_y_file = NULL;

	int width, height, bitdepth;

	ori_yuvfilename = argv[1];
	pre_yuvfilename = argv[2];
	rec_yuvfilename = argv[3];
	width = atoi(argv[4]);
	height = atoi(argv[5]);
	bitdepth = atoi(argv[6]);


	unsigned char* u_buffer = NULL;
	unsigned char* v_buffer = NULL;
	unsigned char* y_buffer = NULL;     //原始图像
	unsigned char* rec_y_buffer = NULL; //重建图像
	unsigned char* pre_y_buffer = NULL; //预测误差

	errno_t err;
	if ((err = fopen_s(&ori_y_file, ori_yuvfilename, "rb")) != 0) {
		cout << "FAIL TO OPEN YUV FILE!";
		exit(1);
	}
	if ((err = fopen_s(&pre_y_file, pre_yuvfilename, "wb")) != 0) {
		cout << "FAIL TO OPEN PRE_YUV FILE!";
		exit(1);
	}
	if ((err = fopen_s(&rec_y_file, rec_yuvfilename, "wb")) != 0) {
		cout << "FAIL TO OPEN REC_YUV FILE!";
		exit(1);
	}

	//开辟空间
	y_buffer = (unsigned char*)malloc(width * height);
	u_buffer = (unsigned char*)malloc(width * height / 4);
	v_buffer = (unsigned char*)malloc(width * height / 4);

	pre_y_buffer = (unsigned char*)malloc(width * height);
	rec_y_buffer = (unsigned char*)malloc(width * height);

	if (y_buffer == NULL || u_buffer == NULL || v_buffer == NULL || pre_y_buffer == NULL|| rec_y_buffer==NULL) {
		cout << "no enought memory\n";
		exit(1);
	}

	//读取原始图像数据
	fread(y_buffer, 1, width * height, ori_y_file);
	fread(u_buffer, 1, width * height / 4, ori_y_file);
	fread(v_buffer, 1, width * height / 4, ori_y_file);

	//得到原始图像的灰度值分布
	FILE* orig;
	
	if ((err = fopen_s(&orig, "ori.txt", "wb")) != 0) {
		cout << "FAIL TO OPEN TXT FILE!";
		exit(1);
	}
	double frequency[256] = { 0 };
	GetFrequency(y_buffer, frequency, height, width);
	for (int i = 0;i < 256;i++)
	{
		fprintf(orig, "%d\t%f\n", i, frequency[i]);
	}
	

	DPCM(y_buffer, pre_y_buffer, rec_y_buffer, width, height, bitdepth);
	PSNR(y_buffer, rec_y_buffer, width, height);

	FILE* pre;

	if ((err = fopen_s(&pre, "pre.txt", "wb")) != 0) {
		cout << "FAIL TO OPEN TXT FILE!";
		exit(1);
	}
	double frequency_pre[256] = { 0 };
	GetFrequency(pre_y_buffer, frequency_pre, height, width);
	for (int i = 0;i < 256;i++)
	{
		fprintf(pre, "%d\t%f\n", i, frequency_pre[i]);
	}
	

	//写入重建图像
	fwrite(rec_y_buffer, width * height, 1, rec_y_file);
	fwrite(u_buffer, width * height / 4, 1, rec_y_file);
	fwrite(v_buffer, width * height / 4, 1, rec_y_file);

	//写入预测图像
	fwrite(pre_y_buffer, width * height, 1, pre_y_file);
	fwrite(u_buffer, width * height / 4, 1, pre_y_file);
	fwrite(v_buffer, width * height / 4, 1, pre_y_file);

	free(y_buffer);
	free(u_buffer);
	free(v_buffer);
	free(pre_y_buffer);
	free(rec_y_buffer);

	fclose(ori_y_file);
	fclose(pre_y_file);
	fclose(rec_y_file);

	return 0;
}