DPCM 压缩系统的实现和分析
一、实验目的
掌握DPCM编解码系统的基本原理。初步掌握实验用C/C++/Python等语言编程实现DPCM
编码器,并分析其压缩效率。
二、实验内容
1.DPCM编解码原理
DPCM是差分预测编码调制的缩写,是比较典型的预测编码系统。在DPCM系统中,
需要注意的是预测器的输入是已经解码以后的样本。之所以不用原始样本来做预测,是
因为在解码端无法得到原始样本,只能得到存在误差的样本。因此,在DPCM编码器中实
际内嵌了一个解码器,如编码器中虚线框中所示。
在一个DPCM系统中,有两个因素需要设计:预测器和量化器。理想情况下,预测器
和量化器应进行联合优化。实际中,采用一种次优的设计方法:分别进行线性预测器和
量化器的优化设计。
2.DPCM编码系统的设计
在本次实验中,我们采用固定预测器和均匀量化器。
量化器采用8比特均匀量化。本实验的目标是验证DPCM编码的编码效率。首先读取一个256级的灰度图像,采用自己设定的预测方法计算预测误差,并对预测误差进行8比特均匀量化
在DPCM编码器实现的过程中可同时输出预测误差图像和重建图像。将预测误差图像写入文件并将该文件输入Huffman编码器,得到输出码流、给出概率分布图并计算压缩比。
将原始图像文件输入Huffman编码器,得到输出码流、给出概率分布图并计算压缩比。最后比较两种系统(1.DPCM+ 熵编码和2. 仅进行熵编码)之间的编码效率(压缩比和图像质量)。压缩质量以PSNR进行计算。
三、实验代码
mian函数:
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include "dpcm.h"
#include <windows.h>
#define u_int8_t unsigned __int8
#define u_int unsigned __int32
#define u_int32_t unsigned __int32
#pragma pack(1)
BITMAPFILEHEADER File_header;
BITMAPINFOHEADER Info_header;
long frameWidth ; /* --width=<uint> */
long frameHeight ; /* --height=<uint> */
unsigned char rectify(unsigned char x)
{
if (x<0) x=0;
if (x>255) x=255;
return x;
}
int BMP2YUV(char *url,char *url_out1,char *url_out2,int num)
{
FILE *bmpFile=NULL;
FILE *pFile=NULL;//存储重建图像值
FILE *eFile=NULL;//存储差值
bool flip = TRUE;/* --flip */
unsigned int i;
unsigned char* bmpBuf = NULL;
unsigned char* pBuf = NULL;//存储重建图像值
//int eBuf;
unsigned char* eBuf_p = NULL;//存储差值
unsigned char* yBuf = NULL;//仅对亮度进行预测编码
unsigned char* uBuf = NULL;
unsigned char* vBuf = NULL;
u_int32_t videoFramesWritten = 0;
bmpFile=fopen(url,"rb+");
pFile=fopen(url_out1,"wb+");
eFile=fopen(url_out2,"wb+");
//判断文件头是否读取正确
if(fread(&File_header,sizeof(BITMAPFILEHEADER),1,bmpFile) != 1)
{
printf("read file header error!");
exit(0);
}
//判断是否为bmp文件
if (File_header.bfType != 0x4D42)
{
printf("Not bmp file!");
exit(0);
}
else
{ printf("this is a bmp\n");
}
//判断信息头是否读取正确
if(fread(&Info_header,sizeof(BITMAPINFOHEADER),1,bmpFile) != 1)
{
printf("read info header error!");
exit(0);
}
// end read header
//判断长宽高是否读取正确,并赋值
if (((Info_header.biWidth/8*Info_header.biBitCount)%4) == 0)
frameWidth= Info_header.biWidth;
else
frameWidth = (Info_header.biWidth*Info_header.biBitCount+31)/32*4;
if ((Info_header.biHeight%2) == 0)
frameHeight = Info_header.biHeight;
else
frameHeight = Info_header.biHeight + 1;
printf("The width is %d\n",frameWidth);
printf("The height is %d\n",frameHeight);
/* open the BMP file */
if (bmpFile == NULL)
{
printf("cannot find bmp file\n");
exit(1);
}
else
{
printf("The input bmp file is open\n"/*, bmpFileName*/);
}
/* open the pfile */
if (pFile == NULL)
{
printf("cannot find p file\n");
exit(1);
}
else
{
printf("The output pfile is open \n"/*, yuvFileName*/);
}
/* open the efile */
if (eFile == NULL)
{
printf("cannot find e file\n");
exit(1);
}
else
{
printf("The output efile is open \n");
}
bmpBuf = (unsigned char*)malloc(frameWidth * frameHeight * 3);
memset(bmpBuf,0,frameHeight*frameWidth*3);
/* get the output buffers for a frame */
yBuf = (unsigned char*)malloc(frameWidth * frameHeight);
//memset(yBuf,0,frameHeight*frameWidth);
pBuf = (unsigned char*)malloc(frameWidth * frameHeight);
eBuf_p =(unsigned char*)malloc(frameWidth * frameHeight);
uBuf = (unsigned char*)malloc((frameWidth * frameHeight)/4 );
vBuf = (unsigned char*)malloc((frameWidth * frameHeight) /4);
/* get an input buffer for a frame */
if (bmpBuf == NULL || yBuf == NULL || pBuf == NULL || eBuf_p == NULL)
{
printf("no enought memory\n");
exit(1);
}
//
void ReadRGB(FILE * pFile,BITMAPFILEHEADER & file_h,BITMAPINFOHEADER & info_h,unsigned char * rgbData);
ReadRGB(bmpFile, File_header, Info_header,bmpBuf);
//while (fread(bmpBuf, 1, frameWidth * frameHeight*3 , bmpFile))
//{
if(RGB2YUV(frameWidth, frameHeight, bmpBuf , yBuf, uBuf, vBuf, flip))
{
printf("error");
return 0;
}
int bits=8;
int scale=512/(1<<bits);
printf("%d\n",scale);
int eBuf ;//存储差值
//一阶
//for( int j=0;j<frameHeight;j++)//j:行
//{
// eBuf=yBuf[j*frameWidth]-128;
// eBuf_p[j*frameWidth]=(eBuf+255)/2;
// yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-128)*2;
// pBuf[frameWidth*j] = (eBuf_p[frameWidth*j] - 128)*2 + 128;
// for(i=1+j*frameWidth;i<frameWidth*(j+1);i++)
// {//i:列
// eBuf=yBuf[i]-yBuf[i-1];
// eBuf_p[i]=(eBuf+255)/2;//预测误差的量化值
// yBuf[i] = (eBuf_p[i] - 128)*2 + yBuf[i-1];
// pBuf[i] = (eBuf_p[i] - 128)*2 + pBuf[i-1];
// }//
//}
//4bit量化
//for( int j=0;j<frameHeight;j++)//j:行
//{
// eBuf=yBuf[j*frameWidth]-128;
// eBuf_p[j*frameWidth]=(eBuf+255)/2;
// yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-127)*2;
// pBuf[frameWidth*j] = (eBuf_p[frameWidth*j] - 127)*2 + 128;
// for(i=1+j*frameWidth;i<frameWidth*(j+1);i++)
// {//i:列
// eBuf=yBuf[i]-yBuf[i-1];
// eBuf_p[i]=(eBuf+255)/2;//预测误差的量化值
// yBuf[i] =(eBuf_p[i] - 127)*2+ yBuf[i-1];
// pBuf[i] = (eBuf_p[i] - 127)*2 + pBuf[i-1];
// }//
//}
//二阶2
for( int j=1;j<frameHeight;j++)//j:行
{
eBuf=yBuf[j*frameWidth]-128;
eBuf_p[j*frameWidth]=(eBuf+255)/2;
yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-255/2)*2;
pBuf[frameWidth*j] = rectify((eBuf_p[frameWidth*j] - 255/2)*2 + 128);
for(i=1;i<frameWidth;i++)
{//i:列
eBuf=yBuf[i+j*frameWidth]-0.5*yBuf[i+j*frameWidth-1]-0.25*yBuf[i+(j-1)*frameWidth]-0.25*yBuf[i+1+(j-1)*frameWidth];
eBuf_p[i+j*frameWidth]=(eBuf+255)/2;//预测误差的量化值
yBuf[i+j*frameWidth] = rectify((eBuf_p[i+j*frameWidth] - 255/2)*2 + 0.5*yBuf[i+j*frameWidth-1]+0.25*yBuf[i+(j-1)*frameWidth]+0.25*yBuf[i+1+(j-1)*frameWidth]);
pBuf[i+j*frameWidth] = rectify((eBuf_p[i+j*frameWidth] - 255/2)*2 + 0.5*pBuf[i+j*frameWidth-1]+0.25*pBuf[i+(j-1)*frameWidth]+0.25*pBuf[i+1+(j-1)*frameWidth]);
}//
}
//一阶0.95
//for( int j=0;j<frameHeight;j++)//j:行
//{
// eBuf=yBuf[j*frameWidth]-128;
// eBuf_p[j*frameWidth]=(eBuf+255)/2;
// yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-255/2)*2;
// pBuf[frameWidth*j] = (eBuf_p[frameWidth*j] - 255 /2)*2 + 128;
// for(i=1+j*frameWidth;i<frameWidth*(j+1);i++)
// {//i:列
// eBuf=yBuf[i]-0.95*yBuf[i-1];
// eBuf_p[i]=(eBuf+255)/2;//预测误差的量化值
// yBuf[i] = (eBuf_p[i] - 127)*2 + 0.95*yBuf[i-1];
// pBuf[i] = (eBuf_p[i] - 127)*2 + 0.95*pBuf[i-1];
// }//
//}
for (i = 0; i < frameWidth*frameHeight; i++)
{
if (pBuf[i] < 16) pBuf[i] = 16;
if (pBuf[i] > 235) pBuf[i] = 235;
}
for (i = 0; i < frameWidth*frameHeight/4; i++)
{
if (uBuf[i] < 16) uBuf[i] = 16;
if (uBuf[i] > 240) uBuf[i] = 240;
if (vBuf[i] < 16) vBuf[i] = 16;
if (vBuf[i] > 240) vBuf[i] = 240;
}
fwrite(pBuf, 1, frameWidth * frameHeight, pFile);
fwrite(uBuf, 1, (frameWidth * frameHeight)/4 , pFile);
fwrite(vBuf, 1, (frameWidth * frameHeight)/4 , pFile);
fwrite(eBuf_p, 1, frameWidth * frameHeight, eFile);
fwrite(uBuf, 1, (frameWidth * frameHeight)/4 , eFile);
fwrite(vBuf, 1, (frameWidth * frameHeight)/4 , eFile);
printf("\r...%d", ++videoFramesWritten);
//}
printf("\n%u %ux%u video frames written\n",
videoFramesWritten, frameWidth, frameHeight);
free(bmpBuf);
free(yBuf);
free(uBuf);
free(vBuf);
/* cleanup */
fclose(bmpFile);
fclose(pFile);
fclose(eFile);
return 0;
}
void ReadRGB(FILE * pFile,BITMAPFILEHEADER & file_h,BITMAPINFOHEADER & info_h,unsigned char * rgbData)
{
bool MakePalette(FILE * pFile,BITMAPFILEHEADER &file_h,BITMAPINFOHEADER & info_h,RGBQUAD *pRGB_out);
unsigned long Loop,iLoop,jLoop,width,height,w,h;
unsigned char mask, *Index_Data,* Data;
if (((info_h.biWidth/8*info_h.biBitCount)%4) == 0)
w = info_h.biWidth;
else
w = (info_h.biWidth*info_h.biBitCount+31)/32*4;
if ((info_h.biHeight%2) == 0)
h = info_h.biHeight;
else
h = info_h.biHeight + 1;
width = w/8*info_h.biBitCount;
height = h;
Index_Data = (unsigned char *)malloc(height*width);
Data = (unsigned char *)malloc(height*width);
fseek(pFile,file_h.bfOffBits,0);
if(fread(Index_Data,height*width,1,pFile) != 1)
{
printf("read file error!\n\n");
exit(0);
}
for ( iLoop = 0;iLoop < height;iLoop ++)
for (jLoop = 0;jLoop < width;jLoop++)
{
Data[iLoop*width+jLoop] = Index_Data[(height-iLoop-1)*width+jLoop];
}
switch(info_h.biBitCount)
{
case 24:
memcpy(rgbData,Data,height*width);
if(Index_Data)
free(Index_Data);
if(Data)
free(Data);
return;
case 16:
if(info_h.biCompression == BI_RGB)
{
for (Loop = 0;Loop < height * width;Loop +=2)
{
*rgbData = (Data[Loop]&0x1F)<<3;
*(rgbData + 1) = ((Data[Loop]&0xE0)>>2) + ((Data[Loop+1]&0x03)<<6);
*(rgbData + 2) = (Data[Loop+1]&0x7C)<<1;
rgbData +=3;
}
}
if(Index_Data)
free(Index_Data);
if(Data)
free(Data);
return;
default:
RGBQUAD *pRGB = (RGBQUAD *)malloc(sizeof(RGBQUAD)*(unsigned long)pow(float(2),info_h.biBitCount));
int temp = sizeof(pRGB);
if(!MakePalette(pFile,file_h,info_h,pRGB))
printf("No palette!\n\n");
for (Loop=0;Loop<height*width;Loop++)
{
switch(info_h.biBitCount)
{
case 1:
mask = 0x80;
break;
case 2:
mask = 0xC0;
break;
case 4:
mask = 0xF0;
break;
case 8:
mask = 0xFF;
}
int shiftCnt = 1;
while (mask)
{
unsigned char index = mask == 0xFF ? Data[Loop] : ((Data[Loop] & mask)>>(8 - shiftCnt * info_h.biBitCount));
* rgbData = pRGB[index].rgbBlue;
* (rgbData+1) = pRGB[index].rgbGreen;
* (rgbData+2) = pRGB[index].rgbRed;
if(info_h.biBitCount == 8)
mask =0;
else
mask >>= info_h.biBitCount;
rgbData+=3;
shiftCnt ++;
}
}
if(Index_Data)
free(Index_Data);
if(Data)
free(Data);
// if(pRGB)
// free(pRGB);
}
}
bool MakePalette(FILE * pFile,BITMAPFILEHEADER &file_h,BITMAPINFOHEADER & info_h,RGBQUAD *pRGB_out)
{
if ((file_h.bfOffBits - sizeof(BITMAPFILEHEADER) - info_h.biSize) == sizeof(RGBQUAD)*pow(float(2),info_h.biBitCount))
{
fseek(pFile,sizeof(BITMAPFILEHEADER)+info_h.biSize,0);
fread(pRGB_out,sizeof(RGBQUAD),(unsigned int)pow(float(2),info_h.biBitCount),pFile);
return true;
}
else
return false;
// free(pRGB);
}
int main(int argc, char **argv)
{
int count=1;
u_int num;
FILE *bmpFile=NULL;
FILE *pFile=NULL;
FILE *eFile=NULL;
char *bmpFileName=NULL;
char *pFileName=NULL;
char *eFileName=NULL;
unsigned char* pBuf = NULL;
unsigned char* eBuf = NULL;
bmpFileName=argv[1];
pFileName=argv[2];
eFileName=argv[3];
num=1;
bmpFile=fopen(bmpFileName,"rb+");
pFile=fopen(pFileName,"wb+");
eFile=fopen(eFileName,"wb+");
BMP2YUV(bmpFileName,pFileName,eFileName,1);
//pBuf = (unsigned char*)malloc((frameWidth * frameHeight)*3/2);
return 0;
}
DPCM.cpp
int RGB2YUV (int x_dim, int y_dim, void *bmp, void *y_out, void *u_out, void *v_out, int flip)
{
static int init_done = 0;
long i, j, size;
unsigned char *r, *g, *b;
unsigned char *y, *u, *v;
unsigned char *pu1, *pu2, *pv1, *pv2, *psu, *psv;
unsigned char *y_buffer, *u_buffer, *v_buffer;
unsigned char *sub_u_buf, *sub_v_buf;
if (init_done == 0)
{
InitLookupTable();//程序一开始,便运行一次亮度计算公式的组件
init_done = 1;
}
// check to see if x_dim and y_dim are divisible by 2,查看是否可以被2整除
if ((x_dim % 2) || (y_dim % 2)) return 1;//如果不能被2 整除,则return1
size = x_dim * y_dim;
// allocate memory
y_buffer = (unsigned char *)y_out;
sub_u_buf = (unsigned char *)u_out;
sub_v_buf = (unsigned char *)v_out;
u_buffer = (unsigned char *)malloc(size * sizeof(unsigned char));//分配指定字节的内存空间
v_buffer = (unsigned char *)malloc(size * sizeof(unsigned char));
if (!(u_buffer && v_buffer))
{
if (u_buffer) free(u_buffer);//如果其中任何一个没有分配内存空间,就释放
if (v_buffer) free(v_buffer);
return 2;
}
b = (unsigned char *)bmp;//b指向bmp文件的首地址的指针
y = y_buffer;
u = u_buffer;
v = v_buffer;
// convert RGB to YUV
if (!flip) {
for (j = 0; j < y_dim; j ++)
{
y = y_buffer + (y_dim - j - 1) * x_dim;
u = u_buffer + (y_dim - j - 1) * x_dim;
v = v_buffer + (y_dim - j - 1) * x_dim;
for (i = 0; i < x_dim; i ++) {
g = b + 1;
r = b + 2;
*y = (unsigned char)( RGBYUV02990[*r] + RGBYUV05870[*g] + RGBYUV01140[*b]);
*u = (unsigned char)(- RGBYUV01684[*r] - RGBYUV03316[*g] + (*b)/2 + 128);
*v = (unsigned char)( (*r)/2 - RGBYUV04187[*g] - RGBYUV00813[*b] + 128);
b += 3;
y ++;
u ++;
v ++;
}
}//倒序
} else {
for (i = 0; i < size; i++)
{
g = b + 1;
r = b + 2;
*y = (unsigned char)( RGBYUV02990[*r] + RGBYUV05870[*g] + RGBYUV01140[*b]);//强制类型转换会产生失真
*u = (unsigned char)(- RGBYUV01684[*r] - RGBYUV03316[*g] + (*b)/2 + 128);//调用函数按照亮度计算公式计算yuv
*v = (unsigned char)( (*r)/2 - RGBYUV04187[*g] - RGBYUV00813[*b] + 128);
b += 3;
y ++;
u ++;
v ++;
}
8bit量化结果:
左侧为原始图像,中间为误差图像,右侧为重建图像。
行huffman编码,输出码流、给出符号概率分布图、计算压缩比如下:
压缩比对比分析
huffman压缩前大小 | huffman压缩后大小 |
---|---|
64 | 58 |
压缩比:1.103
计算PSNR
代码:
def PSNR(pred, gt,height,width):
pred = np.array(pred)
gt = np.array(gt)
pred=pred.reshape(-1,height* width)
gt = gt.reshape(-1, height* width)
res = np.mean((pred - gt) ** 2, axis=1)
res=res.reshape(-1,1)
res=np.sqrt(res)
res[res == 0.0] = 0.01
psnr = 20 * np.log10(255.0 / (res))
# psnr = 20 * np.log10(1.0/ (res))
return psnr
结果:7