图像编码之赫夫曼编码
图像编码
写在前面，这一次的实验相当有难度，完全自己写体会很深刻，把C/C++对于文件读写部分的内容相当全面的覆盖到了，并且也对算法的设计有较高要求。
按照惯例我们还是先说说为什么要对图像进行编码。现在4K逐渐开始普及了，8K的产品也开始销售了，但是这些超高清背后的数据量，你有考虑过嘛？
我们拿现在已经不算稀奇的FullHD（1080p）来做一笔计算：
一幅（帧）图像的字节数：
120分钟的电影：
如果我们用比较好的50G蓝光光盘保存，需要：
张蓝光光碟。
可见，除非Sony宣布蓝光光盘免费，否则看一场电影的门槛可能都要上百了。正因如此，我们才要将图像进行“压缩”，图像能够被压缩的客观条件，就是“冗余”的存在。
一幅图片从不同角度来看会存在多种多样的冗余信息（可以参考这篇科普），我们适当的删减掉不容易被人眼注意到的部分，就可以在不减少信息量的情况下压缩图片的大小。当然这个压缩的具体效果肯定还是与删去的冗余信息量有关的，一幅JPEG压缩的图像在多数情况下都会比原图片看起来“糊”不少，就是因为JPEG的压缩比能够达到26：1左右。
图像编码有非常多的方式，这次的实验内容，就是完成一种基于熵编码的算法——赫夫曼编码。原理这方面不再赘述了，大把的博客以及书籍可供参考，我们直接进入实验内容。
一、实验内容
将一幅给定的灰度图像进行Huffman编解码
众所周知，字数越少，事情越大，这短短一行的实验要求耗费了我三整天的时间。根据我做这个实验过程中遇到的问题，我来将这个实验内容拆分一下：
- 编写函数，获得赫夫曼编码表
- 根据赫夫曼编码，将像素数据用编码表示
- 将编码后的图片数据写入文件
- 读取压缩后的图像文件和赫夫曼编码表
- 利用赫夫曼编码表还原图像文件
二、代码实现与分析
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字节的信息头和文件头，1000字节的调色板，因此像素数据域从第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'，表示的是：
同理，j=12时指向的'0'表示的是：
因此我们就得到了这个通式：
+=
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、利用赫夫曼编码表还原图像文件
这一部分应该还不是最终版，因为我很不满意现在的执行效率，将样例图片还原需要一分钟的时间，这显然是不合格的。
简单说一下目前的思路，我们有了像素数据合并成的一个大字符串data，赫夫曼编码是不定长编码，每次取几个值来核对编码呢？这就需要我们求一下码表中的最短码长和最长码长。
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;
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 = min; k <= max; k++) //从最小码长到最大码长寻找对应码字
{
for (l = 1; l <= 256; l++)
{
if (HCT.length[l] == k) //对码长为k的码字，做比较
{
string temp;
temp += HCT.code[l];
if (!data.substr(index, k).compare(temp)) //注意string的compare方法，相等返回0
{
des.pDataAt(i)[j] = l;
flag = 1; //成功找到解码，标志置1
break;
}
}
}
if (flag == 1) //如果已经找到了像素值，就不再遍历
{
index += k; //下标移动
break;
}
}
if (flag == 0)
{
printf("编码数据有误！\n");
return false;
}
}
}
des.SaveBMPFILE("Recover.bmp");
这个方法简单直观，当然也很慢。显然码长并不是连续的，因此从最短码长遍历到最长码长做了很多无用功，我在想出优化方式后会对博客进行修改。
结语
这一次的实验终于是充满挑战性了，其实这篇博客已经把难点都规避掉了，因为真正的难点在于无从下手，在于走上了错误的解决道路却不自知，我和另外两名同学一起，花了整整三天时间，走遍了岔路才最终得到了结果。当然，收获确实是满满，就是头发又稀了不少。。。
本次实验完整的项目文件与代码可见我的gitee。
