DPCM 压缩系统的实现和分析

读取一个256级的灰度图像，采用左侧/上方预测预测方法计算预测误差，并对预测误差进行各bit量化，DPCM编码器实现的过程中同时输出预测误差图像和重建图像，将该文件Huffman编码器，得到输出码流、给出概率分布图并计算压缩比。将原始图像文件输入Huffman编码器，得到输出码流、给出概率分布图并计算压缩比。最后比较两种系统（1.DPCM+熵编码和2.仅进行熵编码）之间的编码效率（压缩比和图像质量）。压缩质量以PSNR进行计算。

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

2.PSNR
psnr是“Peak Signal to Noise Ratio”的缩写，即峰值信噪比，是一种评价图像的客观标准，它具有局限性，一般是用于最大值信号和背景噪音之间的一个工程项目。单位：dB。

其中，MSE是原图像（语音）与处理图像（语音）之间均方误差。n是每个采样值的比特数。
(以上摘自百度百科）

1.	初始化三个缓存区，分别存放原图像数据original、预测误差predict_error、重建图像数据rebuild；
2.	根据原图像大小分配缓存区，打开原图像、预测误差图像、重建图像；
3.	读取原图像，数据存入缓存区original中；
4.	依次对original中各点灰度值进行统计，得到频次，打开存放原图像灰度频率的文件，逐灰度值计算频率，并写入文件，将得到的频率分布文件绘制成曲线图；
5.	依次取original中灰度数据，判断该处像素是否是图像每行第一个，若是，此处的predict=0,rebuild值等于original中的值；
6.	若不是，此处predict_error灰度值为该原始图像该处灰度值减去重建图像中前一个像素的灰度值，并进行nbit量化；
7.	此处的rebuild灰度值为将predict_error进行nbit反量化后，加上前一像素的rebuild值；
8.	重复5-7，直到取完original中的灰度值；
9.	将predict_error,rebuild中的U、V都赋为128；
10.	将predict_error中的灰度值加上某数，使之全部成为正值；
11.	依次对predict_error中各点灰度值进行统计，得到频次，打开存放原图像灰度频率的文件，逐灰度值计算频率，并写入文件，将得到的频率分布文件绘制成曲线图；
12.	用predict_error,rebuild写预测误差图像，重建图像；
13.	用original,rebuild计算PSNR；
14.	将得到的预测误差图像和原图像一起输入Huffman编码器，得到输出码流，计算压缩比；

2.上方预测
上方预测步骤与左侧预测步骤基本相同，只有5-7步略不一样

5.	依次取original中灰度数据，判断该处像素是否是图像每列第一个，若是，此处的predict=0,rebuild值等于original中的值；
6.	若不是，此处predict_error灰度值为该原始图像该处灰度值减去重建图像中上一个像素的灰度值，并进行nbit量化；
7.	此处的rebuild灰度值为将predict_error进行nbit反量化后，加上上一像素的rebuild值；

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <cmath>
#include <iostream>

#define u_int8_t	unsigned __int8
#define u_int		unsigned __int32
#define u_int32_t	unsigned __int32
#define FALSE		false
#define TRUE		true
using namespace std;

//函数用于写概率分布文件
double getprobability(int frameWidth,int frameHeight,char* FileName,unsigned char* buf) {
	//因UV全为128，所以仅对Y值写出概率分布
	double probailityY[256] = { 0 };
	double HY = 0;
	//打开写入概率分布的文件
	FILE* probaility_File;
	probaility_File = fopen(FileName, "wb");
	int i = 0;
	//依次对各像素上的灰度值进行统计
	for (i = 0; i < frameWidth * frameHeight; i++) {
		int data = buf[i];
		probailityY[data]++;
	}
	//计算频率并写入txt文件
	for (i = 0; i < 256; i++) {
		probailityY[i] = probailityY[i] / (frameWidth * frameHeight);
		//计算熵
		if (probailityY[i] != 0) {
			HY=HY- probailityY[i] * (log(probailityY[i]) / log(2));
		}
		//依次将数据写入txt文件中
		fprintf(probaility_File, "%d\t%f\n", i, probailityY[i]);
	}
	//关闭文件
	fclose(probaility_File);
	return HY;
}

int main(int argc, char** argv)
{
	unsigned int i,j;

	//依次为图像宽高，量化bit数，原始图像文件名，原始图像频率分布文件名，预测误差图像文件名，重建图像文件名，预测误差图像频率分布文件名
	//全部从命令行参数写入
	u_int frameWidth;
	u_int frameHeight;
	u_int quantifybit;
	char* original_FileName = NULL;
	char* original_probability_FileName = NULL;
	char* predict_error_FileName = NULL;
	char* rebuild_FileName = NULL;
	char* predict_error_probability_FileName = NULL;
	frameWidth = atoi(argv[1]);
	frameHeight = atoi(argv[2]);
	quantifybit = atoi(argv[3]);
	original_FileName = argv[4];
	original_probability_FileName = argv[5];
	predict_error_FileName = argv[6];
	rebuild_FileName = argv[7];
	predict_error_probability_FileName = argv[8];

	
	double mse = 0;
	double psnr;
	
	//原始图像文件，预测误差图像文件，重建图像文件
	//频率分布图像文件在getprobability()函数中，所以此处没有定义
	FILE* original_File = NULL;
	FILE* predict_error_File = NULL;
	FILE* rebuild_File = NULL;

	//依次为原始图像数据缓冲区，预测误差图像文件数据缓冲区，重建图像文件缓冲区
	u_int8_t* original_Buf = NULL;
	u_int8_t* predict_error_Buf = NULL;
	u_int8_t* rebuild_Buf = NULL;
	
	
	//打开并判断原始文件是否存在
	original_File = fopen(original_FileName, "rb");
	if (original_File == NULL)
	{
		printf("cannot find yuv file\n");
		exit(1);
	}
	//原始文件存在，打开预测误差图像文件，重建图像文件
	predict_error_File = fopen(predict_error_FileName, "wb");
	rebuild_File = fopen(rebuild_FileName, "wb");
	
	//给3个缓存区分配大小，因为图像为4：2：0，所以大小为宽*高*1.5
	original_Buf = (u_int8_t*)malloc(frameWidth * frameHeight * 1.5);
	predict_error_Buf = (u_int8_t*)malloc(frameWidth * frameHeight * 1.5);
	rebuild_Buf = (u_int8_t*)malloc(frameWidth * frameHeight * 1.5);

	//判断缓存区大小是否分配成功
	if (original_Buf == NULL||predict_error_Buf == NULL || rebuild_Buf == NULL)
	{
		printf("no enought memory\n");
		exit(1);
	}
	//分配成功，读入原始图像文件数据
	fread(original_Buf, 1, frameWidth * frameHeight * 1.5, original_File);
	//调用getprobaility()函数，写出原始图像频率分布文件
	double original_HY=getprobability(frameWidth, frameHeight, original_probability_FileName, original_Buf);
	cout <<"original picture's probability distribution file have written" << endl;
	cout <<"original picture's entropy(Y) is:" << original_HY <<"bit/symbol"<< endl;

	//左向预测，DPCM编码
	for (i = 0; i < frameWidth; i++) {
		for (j = 0; j < frameHeight; j++) {
			if (i == 0) {
				//左向预测，所以每行第一个像素灰度值预测误差为0
				predict_error_Buf[j * frameWidth ] = 0;
				//重建图像中每行第一个像素灰度值为原图像像素值
				rebuild_Buf[j * frameWidth ] = original_Buf[j * frameWidth ];
			}
			else {
				//得到预测误差，并进行量化
				predict_error_Buf[j * frameWidth + i] = (original_Buf[j * frameWidth + i] - rebuild_Buf[j * frameWidth + i - 1])/(256 / pow(2, quantifybit - 1));
				//对预测误差进行反量化并得到重建图像像素灰度值
				rebuild_Buf[j * frameWidth + i] = rebuild_Buf[j * frameWidth + i - 1] + predict_error_Buf[j * frameWidth + i]*(256 / pow(2, quantifybit - 1));
			}
		}

	}
	
	//写出预测误差图像
	//预测误差有负值，这里将其全部变为正值
	for (i = 0; i < frameWidth * frameHeight; i++) {
		predict_error_Buf[i] = predict_error_Buf[i] + pow(2,quantifybit-1);
	}
	//黑白图像，UV全部赋为128
	for (i = frameWidth * frameHeight; i < frameWidth * frameHeight * 1.5; i++) {
		predict_error_Buf[i] = 128;
	}
	fwrite(predict_error_Buf, 1, (frameWidth * frameHeight) * 1.5, predict_error_File);
	cout << quantifybit << "bit quantified predict_error_picture have written" << endl;
	//调用getprobaility()函数，写出预测误差图像频率分布文件
	double predict_error_HY=getprobability(frameWidth, frameHeight, predict_error_probability_FileName, predict_error_Buf);
	cout<< quantifybit << "bit quantified predict_error_picture's probability distribution file have written" << endl;
	cout << quantifybit << "bit predict_error_picture's entropy(Y) is:" << predict_error_HY << "bit/symbol" << endl;

	//写出重建图像文件
	//黑白图像，UV全部赋为128图像
	for (i = frameWidth * frameHeight; i < frameWidth * frameHeight * 1.5; i++) {
		rebuild_Buf[i] = 128;
	}
	fwrite(rebuild_Buf, 1, (frameWidth * frameHeight) * 1.5, rebuild_File);
	cout << quantifybit << "bit quantified rebuild_picture have written" << endl;

	//计算图像压缩质量（PSNR）
	for (i = 0; i < frameWidth * frameHeight * 1.5; i++) {
		mse = mse + pow(original_Buf[i] - rebuild_Buf[i], 2);
	}
	mse = mse / (frameWidth * frameHeight*1.5);
	double a = pow(pow(2, quantifybit)-1,2);
	psnr =10 * log10(a/ mse);

	cout << quantifybit<<"bit quantifying encoded and the PSNR is:"<<psnr <<"dB"<<endl;

	//关闭文件，释放缓存空间
	free(original_Buf);
	free(predict_error_Buf);
	free(rebuild_Buf);
	fclose(original_File);
	fclose(predict_error_File);
	fclose(rebuild_File);
	return(0);
}

2.上方预测（代码基本与左侧预测相同，仅DPCM编码部分不同，所以仅这一部分）

for (i = 0; i < frameWidth; i++) {
		//上方预测，所以每列第一个像素灰度值预测误差为0
		predict_error_Buf[i] = 0;
		//重建图像中每列第一个像素灰度值为原图像像素值
		rebuild_Buf[i] = original_Buf[i];
		for (j = 1; j < frameHeight; j++) {
			//得到预测误差，并进行量化
			predict_error_Buf[j * frameWidth + i] = (original_Buf[j * frameWidth + i] - rebuild_Buf[(j - 1) * frameWidth + i]) / (256 / pow(2, quantifybit - 1));
			//对预测误差进行反量化并得到重建图像像素灰度值
			rebuild_Buf[j * frameWidth + i] = rebuild_Buf[(j-1)* frameWidth + i] + predict_error_Buf[j * frameWidth + i] * (256 / pow(2, quantifybit - 1));
		}
	}

3.用python绘制频率分布曲线图（以原始图像的绘制为例）

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate

def loadData(flieName):
    a = np.loadtxt(flieName)
    X = a[:, 0]
    y = a[:, 1]
    return X, y  # X,y组成一个元组，这样可以通过函数一次性返回


if __name__ == '__main__':
    X, y = loadData('original_probability.txt')
    plt.plot(X, y)
    plt.xlabel('Y')
    plt.ylabel('probability')
    plt.title('original_picture')
    # 保存图片
    plt.savefig('original_probability.jpg')
    plt.show()

4.Huffman编码器命令行

1.原图像

2.8bit
左侧预测：
命令行参数：

256 256 8 Lena256B.yuv original_probability.txt 8bit_predict_error.yuv 8bit_rebuild.yuv 8bit_predict_error_probaility.txt

（从左至右依次是原图像，预测误差图像，重建图像，下同）

上方预测：
命令行参数：

256 256 8 Lena256B.yuv original_probability.txt 8bit_uppredict_error.yuv 8bit_uprebuild.yuv 8bit_uppredict_error_probaility.txt

3.4bit
左侧预测：
命令行参数：

256 256 4 Lena256B.yuv original_probability.txt 4bit_predict_error.yuv 4bit_rebuild.yuv 4bit_predict_error_probaility.txt

上方预测：
命令行参数：

256 256 4 Lena256B.yuv original_probability.txt 4bit_uppredict_error.yuv 4bit_uprebuild.yuv 4bit_uppredict_error_probaility.txt

4.2bit
左侧预测：
命令行参数：

256 256 2 Lena256B.yuv original_probability.txt 2bit_predict_error.yuv 2bit_rebuild.yuv 2bit_predict_error_probaility.txt

上方预测：
命令行参数：

256 256 2 Lena256B.yuv original_probability.txt 2bit_uppredict_error.yuv 2bit_uprebuild.yuv 2bit_uppredict_error_probaility.txt

5.1bit
左侧预测：
命令行参数：

256 256 1 Lena256B.yuv original_probability.txt 1bit_predict_error.yuv 1bit_rebuild.yuv 1bit_predict_error_probaility.txt

上方预测：
命令行参数：

256 256 1 Lena256B.yuv original_probability.txt 1bit_uppredict_error.yuv 1bit_uprebuild.yuv 1bit_uppredict_error_probaility.txt

1.图像质量（PSNR）

可以看出8bit量化得出的重建图像基本与原图像相同，质量远大于其他bit量化得到的重建图像。

2.压缩比

可以看到对于这副图像，DPCM+熵编码（nbit量化后再进行哈夫曼编码）的压缩比大于直接进行熵编码（原图像直接进行哈夫曼编码）。