图像编码之赫夫曼编码

• 图像编码
• 一、实验内容
• 二、代码实现与分析
• 1、编写函数，获得赫夫曼编码表
• 2、根据赫夫曼编码，将像素数据用编码表示
• 3、将编码后的图片数据写入文件
• 4、读取压缩后的图像文件和赫夫曼编码表
• 4.1 读取赫夫曼编码表
• 4.2 读取图像数据
• 5、利用赫夫曼编码表还原图像文件
• 结语

图像编码

写在前面，这一次的实验相当有难度，完全自己写体会很深刻，把C/C++对于文件读写部分的内容相当全面的覆盖到了，并且也对算法的设计有较高要求。

按照惯例我们还是先说说为什么要对图像进行编码。现在4K逐渐开始普及了，8K的产品也开始销售了，但是这些超高清背后的数据量，你有考虑过嘛？
我们拿现在已经不算稀奇的FullHD（1080p）来做一笔计算：
一幅（帧）图像的字节数： 1920 ∗ 1080 ∗ 3 = 6.22 M 1920*1080*3=6.22M 1920∗1080∗3=6.22M
120分钟的电影： 6.22 M ∗ 60 ∗ 60 ∗ 120 = 2687 G 6.22M*60*60*120=2687G 6.22M∗60∗60∗120=2687G
如果我们用比较好的50G蓝光光盘保存，需要： 2687 / 50 = 54 2687/50=54 2687/50=54张蓝光光碟。
可见，除非Sony宣布蓝光光盘免费，否则看一场电影的门槛可能都要上百了。正因如此，我们才要将图像进行“压缩”，图像能够被压缩的客观条件，就是“冗余”的存在。
一幅图片从不同角度来看会存在多种多样的冗余信息（可以参考这篇科普），我们适当的删减掉不容易被人眼注意到的部分，就可以在不减少信息量的情况下压缩图片的大小。当然这个压缩的具体效果肯定还是与删去的冗余信息量有关的，一幅JPEG压缩的图像在多数情况下都会比原图片看起来“糊”不少，就是因为JPEG的压缩比能够达到26：1左右。
图像编码有非常多的方式，这次的实验内容，就是完成一种基于熵编码的算法——赫夫曼编码。原理这方面不再赘述了，大把的博客以及书籍可供参考，我们直接进入实验内容。

一、实验内容

将一幅给定的灰度图像进行Huffman编解码

众所周知，字数越少，事情越大，这短短一行的实验要求耗费了我三整天的时间。根据我做这个实验过程中遇到的问题，我来将这个实验内容拆分一下：

1. 编写函数，获得赫夫曼编码表
2. 根据赫夫曼编码，将像素数据用编码表示
3. 将编码后的图片数据写入文件
4. 读取压缩后的图像文件和赫夫曼编码表
5. 利用赫夫曼编码表还原图像文件

二、代码实现与分析

1、编写函数，获得赫夫曼编码表

别看这是唯一一个和标题名字重合的步骤，实际上却是最简单的，当然这也是因为我上学期学习数据结构时就已经完成了赫夫曼树、赫夫曼编码的函数，这里我只进行了简单的修改，便得到了码表。
简单的讲一下代码，为什么这么写可以参考我老师的mooc

typedef struct {
	float weight;	//结点权值
	int parent, lchild, rchild;	//指向双亲和左右孩子的指针
}HTNode,*HuffmanTree;	//动态分配数组赫夫曼树

typedef char **HuffmanCode;	//动态分配数组存储赫夫曼编码

typedef struct {
	HuffmanCode code;	//码字数组，从下标1开始
	int *length;		//码长数组，从下标1开始
}HuffmanCodeTable;	//赫夫曼编码表

上面是结构体的定义，我们的思路是先构造赫夫曼树，再根据赫夫曼树从根节点开始走到各个结点的路径得到赫夫曼编码(左0右1)。上这门课的时候要求使用C实现，实际上用C++的string类型可能更加方便。

void HuffmanCoding(HuffmanTree &HT, HuffmanCode &HC, HuffmanCodeTable &HCT, float *w, int n)
{	//w存放n个字符的权值，构造赫夫曼树HT，并求出w中n个数据对应的赫夫曼编码表HCT
	int i, j, m, s1, s2;
	if (n <= 1)
		return;
	HTNode *p;	//指针
	m = 2 * n - 1;	//m是总的结点数
	HT = (HuffmanTree)malloc((m + 1) * sizeof(HTNode));  //0号单元作为标记
	for (p = HT + 1, i = 1; i <= n; i++, p++, w++)
	{
		*p = { *w,0,0,0 };	//叶子结点初始化
	}
	for (i = n + 1; i <= m; i++, p++)
	{
		*p = { 0,0,0,0 };	//非叶子结点初始化
	}
	for (i = n + 1; i <= m; i++)
	{
		Select(HT, i - 1, s1, s2);	//寻找最小值下标s1和次小值下标s2
		HT[s1].parent = i;
		HT[s2].parent = i;
		HT[i].lchild = s1;
		HT[i].rchild = s2;
		HT[i].weight = HT[s1].weight + HT[s2].weight;
	}	//赫夫曼树建立完毕

	//求赫夫曼编码的准备工作
	HC = (HuffmanCode)malloc((n+1) * sizeof(char*));  //分配n个字符串编码的头指针向量
	char *cd = (char*)malloc(n * sizeof(char)); //分配求当前编码的工作空间
	cd[n - 1] = '\0';	//从右向左逐位存放编码，首先存放编码结束符

	int start, c, f;
	for (i = 1; i <= n; i++)	//按照w中的顺序给出赫夫曼编码
	{
		start = n - 1;	//初始化编码起始指针
		for (c = i, f = HT[i].parent; f != 0; c = f, f = HT[f].parent)	 //从叶子结点到根结点逆向求编码
		{
			if (HT[f].lchild == c)
				cd[--start] = '0';	//左分支标0
			else
				cd[--start] = '1';	//右分支标1
		}
		HC[i] = (char*)malloc((n - start) * sizeof(char));	//为第i个编码分配空间，(n-start)为编码长度
		strcpy(HC[i], &cd[start]);	//从cd复制编码串到HC
	}
	free(cd);

	//对赫夫曼编码的码字求长度，得到赫夫曼编码表
	HCT.code = HC;
	HCT.length = (int*)malloc(sizeof(int)*(n+1));
	int count;
	for (i = 1; i <= n; i++)
	{
		count = 0;	//实际上可以用strlen()
		for (j = 0; HC[i][j] != '\0'; j++)
		{
			count += 1;
		}
		HCT.length[i] = count;
	}
}

赫夫曼编码中有一个Select函数是用来求所有没有父结点的结点中权值最小的两个，从而构成新的结点。我的实现写复杂了，想参考可以见文末传送门，建议写一个排序来实现，效率会更高。

2、根据赫夫曼编码，将像素数据用编码表示

上面的函数只是能够实现赫夫曼编码，实际上我们还没有开始调用，下面来说一下对图像调用赫夫曼编码的想法。我们压缩压的是像素数据域的内容，显然，如果我们把出现次数多的像素值赋予短编码，出现次数少的像素值赋予长编码，那么从均值的角度看，我们还是成功的减小了总长度，也就实现了压缩的目的。
因此，我们应当统计图像中各个像素值（我们这里只考虑灰度图像，故为灰度值）出现的频率，然后对频率进行赫夫曼编码。

	for (i = 0; i < src.imagew; i++)	//统计像素出现次数
		for (j = 0; j < src.imageh; j++)
		{
			count_pixel[src.pDataAt(i)[j]] += 1;
		}
	
	for (i = 0; i < 256; i++)	//计算像素出现频率
	{
		p[i] = (float)count_pixel[i] / (float)n;
	}

	HuffmanCoding(HT, HC, HCT, p, 256);	//赫夫曼编码函数

这样我们就成功的获得了HCT，赫夫曼编码表，它由一个码字数组和一个码长数组构成（均从下标1开始），将其写入文件后，如下：（使用'\t'分隔码字与码长）

下面我们用它来将像素数据域的内容编码。不过在转换编码前，我们应该先做一个思考，应当用什么格式来存储编码后的内容呢？
首先我们要明确，写入文件一定是以二进制形式完成的，因为ASCII码表示的'0'和'1'本身就占一个字节，如果以文本格式输出则完全是与压缩的初衷背向而驰。所以，我们要选一个很容易转换为二进制的格式，字符数组便成了首选。
其次，char*也可以表示字符数组，string也可以表示字符数组，具体用哪个更好呢。我们需要将所有像素编码后的结果合并到一个大字符串中，因此对于char*而言，我们要用strcat(des,src)，对于string而言，我们可以用s1+=s2。经过我的实测，strcat执行时间超过40s，string的+=只需要2s，这么大的差异也确实出乎我的意料，不过这篇博客中的测试也显示了strcat的效率不尽人意。

string data;
for (i = 0; i < src.imagew; i++)
	for (j = 0; j < src.imageh; j++)
	{
		data += HCT.code[src.pDataAt(i)[j] + 1];
	}

3、将编码后的图片数据写入文件

好了，编码已经完成了，下面需要把压缩后的图片数据写入文件了。首先我们把信息头、文件头、调色板的数据先写入：

//下面开始将信息头、文件头，这部分没有编码
	FILE *f;
	f = fopen("encode_Image.bmp","w+b");
	if (f == NULL)
		return false;
	BITMAPFILEHEADER fh;
	BITMAPINFOHEADER ih;
	memset(&ih, 0, sizeof(BITMAPINFOHEADER));
	fh.bfType = 0x4d42;
	fh.bfReserved1 = fh.bfReserved2 = 0;
	fh.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + ((src.iYRGBnum == 1) ? 256 * sizeof(RGBQUAD) : 0);
	ih.biSize = 40;
	ih.biPlanes = 1;
	ih.biWidth = src.imagew;
	ih.biHeight = src.imageh;
	ih.biBitCount = 8 * src.iYRGBnum;
	int w4b = (src.imagew * src.iYRGBnum + 3) / 4 * 4;
	ih.biSizeImage = ih.biHeight * w4b;
	fh.bfSize = fh.bfOffBits + ih.biSizeImage;
	fwrite(&fh, sizeof(BITMAPFILEHEADER), 1, f);
	fwrite(&ih, sizeof(BITMAPINFOHEADER), 1, f);
	if (src.iYRGBnum == 1)	//因为我们只处理灰度图像，所以这里会需要写入调色板
		fwrite(src.palette, sizeof(RGBQUAD), 256, f);
	fclose(f);

下面，我们要处理像素数据域的数据了。刚才第二步中我们把用赫夫曼编码表示的像素数据存到了string data;这个大字符串中，下面我们要把这个’1’、'0’构成的字符串，当成10二进制码写入文件了。
要完成这个操作还是有些困难的，我在尝试将字符转换为bitset无果后，还是去搜寻了一番大神的博客，真的被我找到了一个可以实现我们目标的函数，这里放上博客地址，感谢这位大大。
需要注意的是，我们是在写完信息头、文件头等信息的文件中续写像素数据，因此需要将"wb"模式改为"ab"模式。

writeBit(data, "encode_Image.bmp");	//记得修改writeBit参数为string

可以用sublime text等软件查看其数据：（灰度图像有54字节的信息头和文件头，1024字节的调色板，因此像素数据域从第45行倒数第二个字节开始，如果你和我一样用的1.bmp这张图，其值应为4c）

至此，我们已经完成了对原始图像的赫夫曼编码。

4、读取压缩后的图像文件和赫夫曼编码表

之后是解码的操作，首先我们要从文件读回赫夫曼编码表，以及图像数据，再根据编码表还原出原始数据。

4.1 读取赫夫曼编码表

回顾一下我们写的编码表文件，其格式为一行之中，码字'\t'码长'\n'，我们需要读取这个文件来恢复出一个HuffmanCodeTable HCT;
根据上述的需求，我们可以按行读取，将遇到'\t'之前的部分赋值给码字，将'\t'之后的部分，转换为数字，赋值给码长。代码如下：

bool readHuffmanTable(const char* cFilename,HuffmanCodeTable &HCT)	//从文件读取赫夫曼码表
{
	HCT.code = (HuffmanCode)malloc((256 + 1) * sizeof(char*));	//为HCT分配空间
	HCT.length = (int*)malloc((256 + 1) * sizeof(int));
	memset(HCT.length, 0, 257 * sizeof(int));

	FILE *f;
	f = fopen(cFilename, "r");
	if (f == NULL)
		return false;

	char buffer[50];
	int i=1, j=0, length;
	while (fgets(buffer, 50, f) != NULL)
	{
		length = strlen(buffer);
		HCT.code[i] = (char*)malloc(sizeof(char)*(length+1));
		for (j = 0; j < length; j++)
		{
			if (buffer[j] == '\t')
			{
				HCT.code[i][j] = '\0';
				while (buffer[++j] != '\n')
				{
					HCT.length[i] += (int)((buffer[j] - '0')*(pow((double)10, (double)(length - (j + 1) - 1))));
				}
				break;
			}
			else
				HCT.code[i][j] = buffer[j];
		}
		i++;
	}
	return true;
}

可以注意到，对于码长的赋值，我自己推导了一个公式。举一个简单的例子，很容易搞清楚这个公式。
假设，某一行内容如下：1111100110 10，那么该行的长度length为14（'\t'和结尾的'\n'不要忘记）。读到'\t'时的下标j为10，说明下一位开始表示的就是码长了。但是码长不一定只有一位，因此我们写一个内部循环来求码长。

while (buffer[++j] != '\n')

j=11时指向的是1，这个1应该是个位、十位还是百位呢（如果有的话，实际上百位不太可能）？显然，我们可以通过length与j的值来推导出指数。j=11时指向的'1'，表示的是：
1 0 ( 14 − ( 11 + 1 ) − 1 ) = 1 0 1 10^{(14-(11+1)-1)}=10^{1} 10(14−(11+1)−1)=101
同理，j=12时指向的'0'表示的是：
1 0 ( 14 − ( 12 + 1 ) − 1 ) = 1 0 0 10^{(14-(12+1)-1)}=10^{0} 10(14−(12+1)−1)=100
因此我们就得到了这个通式：
H C T . l e n g t h [ i ] HCT.length[i] HCT.length[i] += 1 0 ( l e n g t h − ( j + 1 ) − 1 ) 10^{(length-(j+1)-1)} 10(length−(j+1)−1)

4.2 读取图像数据

就像写入时一样，我们将信息头、文件头、调色板先读取出来，确定图像的长、宽、位数等基本信息。

	FILE *f;
	if ((f = fopen(cFilename, "r+b")) == NULL)
	{
		return false;
	}
	//文件头、信息头的读取，和LoadBMPFILE函数一致
	BITMAPFILEHEADER fh;	//文件头
	BITMAPINFOHEADER ih;	//信息头
	fread(&fh, sizeof(BITMAPFILEHEADER), 1, f);		//从文件中读取文件头
	if (fh.bfType != 0x4d42)	//如果不是"BM"文件，直接退出
	{
		fclose(f);
		return false;
	}
	fread(&ih, sizeof(BITMAPINFOHEADER), 1, f);		//从文件中读取信息头
	if ((ih.biBitCount != 8) && (ih.biBitCount != 24))	//如果不是8位灰度图，或24位彩图，则退出
	{
		fclose(f);
		return false;
	}
	des.iYRGBnum = ih.biBitCount / 8;	//1为灰度图，3为彩图
	des.imagew = ih.biWidth;
	des.imageh = ih.biHeight;
	if (!des.AllocateMem())	//如果分配内存失败，退出
	{
		fclose(f);
		return false;
	}
	if (des.iYRGBnum == 1)	//灰度图需要调色板
		fread(des.palette, sizeof(RGBQUAD), 256, f);
	//信息头、文件头读取完毕

下面我们要处理像素数据了，之前写入时我们把字符串以二进制输出，现在我们要把这些二进制的01重新存为字符串。好在这学期的另一门课《计算机组成原理》中我曾经做过将float型按位输出以检验是否符合IEEE754标准的实验，因此我有一个函数可以将输入的东西以二进制输出。

char* Byte2Binary(unsigned char buffer)	//将一个字节转换成8bit保存到字符串
{
	char *s = (char*)malloc(sizeof(char) * 9);	//返回字符串，记得给'\0'预留位置
	s[0] = '\0';
	unsigned char *p, ch;
	int i, j;
	p = (unsigned char *)(&buffer);
	for (i = sizeof(buffer) - 1; i >= 0; i--)
	{
		ch = *(p + i);
		for (j = 0; j < 8; j++)
		{
			if ((ch << j) & 0x80)	//0x80即128，二进制表示1000 0000，这里操作目的是取最高位
				s[j] = '1';
			else
				s[j] = '0';
		}
		s[j] = '\0';
	}
	return s;
}

用这个函数，我们只需要按一个字节的大小(恰好是一个unsigned char)来读取文件，再将其返回的字符数组拼接起来即可。

	unsigned char buf;	//注意char是1字节，char*是4字节
	int rc;
	string data;	//获得编码数据
	while ((rc = fread(&buf, sizeof(unsigned char), 1, f) != 0))
	{
		data += Byte2Binary(buf);	//按字节读入，按位输出，存入字符串data
	}

5、利用赫夫曼编码表还原图像文件

如果看过我最早版本的代码，并且运行过，那么你一定还记得这里我解码一幅57kb的图像耗时都高达40s。这里我告诉大家问题所在：字符串的相关操作！
最初我为了偷懒，用了string的自带函数substr(index,length)，然后与码表中的码字compare()。这种算法的最大问题就在于substr需要不断的开辟空间（说不断其实是不准确的，具体的解释可以自行查找），因此不仅占用空间，而且相当费时。下面的版本中，我老老实实的逐字符比较，并且极力避免有关出现string类型的构造，最终达到了5s出结果的成绩。代码如下：

bool decode(HXLBMPFILE &des,HuffmanCodeTable &HCT,string data)	//解码函数
{
	int i, j, k;
	int min = HCT.length[1], max = HCT.length[1];	//求最短码长和最大码长
	for (i = 1; i <= 256; i++)
	{
		if (HCT.length[i] < min)
		{
			min = HCT.length[i];
		}
		if (HCT.length[i] > max)
		{
			max = HCT.length[i];
		}
	}

	int index = 0;	//字符串的下标标识
	int flag;	//匹配标志
	int data_length = data.length();	//数据总长度
	for (i = 0; i < des.imagew; i++)	// 解码，给每个像素赋值
		for (j = 0; j < des.imageh; j++)
		{
			flag = 0;	//标志置0
			if ((data_length - (index + 1)) >= min)	//当剩余字符数大于最小码长时
			{
				for (k = 1; k <= 256; k++)	//遍历码表
				{
					int temp = index;	//备份初始下标
					for (int n = 0; n < HCT.length[k]; n++)
					{
						if (data[index] == HCT.code[k][n])	//如果匹配，下标移动，标志置1
						{
							flag = 1;
							index++;
						}
						else
						{
							index = temp;	//不匹配则下标回到初始位置
							flag = 0;	//标志回置0
							break;
						}
					}
					if (flag == 1)
					{
						des.pDataAt(i)[j] = k;
						break;
					}
					
				}
				if (flag == 0)
				{
					printf("编码数据与码表无法匹配！\n");
					return false;
				}
			}
		}
	return true;
}

结语

这一次的实验终于是充满挑战性了，其实这篇博客已经把难点都规避掉了，因为真正的难点在于无从下手，在于走上了错误的解决道路却不自知，我和另外两名同学一起，花了整整三天时间，走遍了岔路才最终得到了结果。当然，收获确实是满满，就是头发又稀了不少。。。
本次实验完整的项目文件与代码可见我的gitee。