问题：额外都到的字母是什么？画出接收到每一个字母之后的自适应哈夫曼树
比较JPEG与GIF压缩之间的不同，并回答对某一图片进行压缩哪一种效果更好，并验证你的答案
实验目的
- 比较自适应哈夫曼编码与哈夫曼编码的区别
- 熟识JPEG压缩的流程以及GIF的LZW编码
实验内容
第一题：
问题1：自适应哈夫曼编码与原始的哈夫曼编码相比有什么优点？
- 哈夫曼编码需要有有关信息源的先验统计知识，也就是需要提前知道信息源的所有编码信息，是一种全局式的编码算法，但是实际上这样的信息通常是很难获得的，因为数据到达之前都是未知的。
- 同时原始的哈夫曼编码是一种静态的算法，一个哈夫曼表只适用于某一状态的的数据，当数据出现一点点变化的时候就会出现不适用的现象。
- 自适应哈夫曼压缩算法的统计数字是随着数据流的到达而动态地收集和更新的，概率不再是基于先验知识而是基于到目前为止实际都到的数据。它允许在符号正在传输时构建代码，允许一次编码并适应数据中变化的条件，即随着数据流的到达，动态地收集和更新符号的概率（频率）。一遍扫描的好处是使得源程序可以实时编码，但由于单个丢失会损坏整个代码，因此它对传输错误更加敏感。
- 在哈夫曼编码中，有个缺点是除了压缩后的资料外，它还得传送机率表给解码端，否则解码端无法正确地做解码的工作。如果想要压缩好一点，必须有更多的统计资料，但同时必须要送出更多的统计资料到解压缩端。而适应性编码可以利用已经读过的资料机动的调整哈夫曼树。适应性哈夫曼编码中，算法FGK的基本原则是兄弟性质
问题2：初始的编码是a = 00，b = 01，c = 10，d = 11，一个符号的第一次发送，则必须首先发送特殊符号NEW。如图所示是发送字母aabb之后的自适应哈夫曼编码。之后，解码器接收到的后续的几个字母的编码位流是01010010101
问题：额外都到的字母是什么？画出接收到每一个字母之后的自适应哈夫曼树
由于上面两道题目是相互关联的，因此可以合并为一题进行解答：
1、 可知当前的a、b、new的编码分别为1，01，00，所以可知解码器首先接收到的编码为01，对应的是b，所以自适应哈夫曼树变化为：
2、B的权重增加1，所以b的父节点的权重变为3，从而大于2，因此会与a节点进行交换，也就是发生以下的变换
3、因此下一个编码01也就是a，所以哈夫曼树中a的权重增加1，如下
4、下一个编码为00，也就是NEW的编码，所以接下来的编码是新的字母，10，也就是c的初始编码，添加到哈夫曼树中会发生如下的变化：
5、可以看到a的父节点的权重大于b的权重，因此就会发生交换，从而如下：
6、因此下一个编码是101，也就是c，所以哈夫曼树变为
因此解码器接收到的字母为：bacc，自适应哈夫曼树的变化如上面的步骤显示：
第二题：
比较JPEG与GIF压缩之间的不同，并回答对某一图片进行压缩哪一种效果更好，并验证你的答案
结论：
- 选择JPEG
- 一般情况下，JPEG的压缩率比GIF小。JPEG的压缩是一个范围，因为可以通过给量化表前面添加一个值来改变压缩率，因此JPEG的压缩率可以比GIF的大，也可以比GIF的小。
- 一般情况下JPEG的视觉效果会比GIF的好。
理论原因如下：
- JPEG图片是以24为颜色存储单个光栅图像。是与平台无关的格式。支持最高级别的压缩。可以通过改变量化表的参数提高与降低JPEG文件压缩的级别。JPEG压缩比可以高达100：1，而且JPEG格式在10：1到20：1的比例下可以轻松压缩文件并且图片的质量仍可以保持良好。JPEG压缩可以很好的处理写实摄影作品。但是对于颜色比较少而且对比度很强烈的，如卡通图片、实心边框或纯色区域大的简单作品，JPEG压缩无法提供理想的结果。优势压缩比会降到5：1，并且严重损失了图片的完整性。这一种损失的原因是JPEG压缩方案可以很好的压缩类似的色调。但是不能够很好地处理亮度的强烈差异或处理纯色区域。
- JPEG压缩在从颜色空间转化到二次采样的时候压缩比就已经达到了2：1，而在dct与量化变化之后数据就已经集中于左上角，形成的矩阵中含有很多的0，这时候采用游长编码可以降低很多的数据量，同时因为整个图片的dc一般都不会相差很大，所以使用DPCM进行压缩也降低了很多的数据量。最后在前面RLC与DPCM的基础上进行哈夫曼编码使得数据量再次大幅度下降。另外哈夫曼编码是一种使用变长编码表对源符号进行编码的，这也就是意味着压缩后的数据还需要添加变长编码表，所以这是哈夫曼编码所带的不足之处。
- JPEG的数据的损失主要是在dct变换与量化的步骤中，这是因为在量化的步骤中是进行取整的，所以很多的相似的数据经过量化之后会变成一样的数据，发生归一，在经过逆变换之后就相同了，也就是发生了损失，另外dct变换中只保留了低阶的数据，较高阶的数据则被丢弃。
- GIF图片文件的数据是经过压缩处理过的，也就是图片中的RGB值只有256种，相当于在一种256色图片基础上的索引类型的的图片。其数据存放在LZW算法的压缩编码流中，文件数据格式包括文件头以及文件数据，文件头包括GIF签名，屏幕描述符、全局色图等。在存储的时候使用一种公共的索引表，并将图片中使用的颜色提取出来，组成一个调色盘，而在存储图片点阵的时候只需要存储每一个点在调色盘中的索引值。将调色盘存储在文件头之中，作为共用的信息。这样的话一个256色的调色盘，24bit的颜色就可以使用8bit来表达了。这样的压缩率为3：1了。另外因为是256色图片的基础上，所以对于颜色总的数目比较少的图片压缩效果会比较好。
- LZW编码是自适应码，不需要知道信源符号的概率分布。同时是字典码，在构造的过程中根据信源序列中出现的信源字符串构造一个字符串表，用定长的单词序号作为该单词的码字，随着码字的字长增加编码效率趋近于1.另外压缩的速率比较快。但是相比于上面的混合压缩编码方式LZW的压缩效果并不能占优势。
程序实现代码如下：
- 图片的扩充：
因为后面图片需要分割为8×8的小块作为dct变换处理的矩阵，所以图片的像素长与高至少需要扩充为8的倍数。另外因为后面的颜色空间变换中U、V与Y的数目比值为1：2，因此需要改为至少为16的倍数。扩充的像素点的值以0代替。因此实现的代码如下：
# 因为UV是二次采样，而且后面需要dct变换，所以需要除以16
# 输入：原图片
# 输出：扩展后的图片（长宽均可以被16整除）
def spread_image(image):
# 计算需要扩展的col的数目
col = 16 - width % 16
# 计算需要扩展的rwo的数目
row = 16 - height % 16
# 生成一张新的图片，宽度加上原本的宽度加上需要扩展的宽度
spread_Image = Image.new('RGB', (width+col, height+row))
new_pixel = spread_Image.putpixel
# 获取图片中的像素点的数据复制（没有进行复制的则默认为0），也就是以0扩展
old_pixel = image.load()
# 进行每个像素上值的
for i in range (width):
for j in range(height):
new_pixel((i,j), old_pixel[i,j])
return spread_Image- 颜色空间的变换：
JPEG采用的是YUV420的颜色空间，而BMP采用的是RGB的颜色空间，因此在进行压缩之前需要先进行颜色的转换。其中YUV包括YCrCb等多种格式，其中Y代表的是亮度，U、V代表的是色度。JPEG常用的RGB到YUV的转换公式如下：
Y = 0.229R + 0.587G + 0.114B
U = -0.169R - 0.331G + 0.5B + 128
V = 0.5R - 0.419G - 0.081B + 128
运算的过程可以使用矩阵实现。
另外还需要对YUV颜色空间进行二次采样，使用的是YUV420，也就是将图片逻辑上划分为2×2的方格（因为前面进行了图片的扩充，因此不存在不能够被整除的情况），每一个方格保留4个亮度值，然后保留2个色度值，这两个色度值可以是两个U或者两个V或者两者各一个，这与YUV4：1：1的区别是YUV4：1：1是一个1×4的方格内的取样。所以颜色空间转换以及采样的代码实现如下：
# 将RGB转化为YUV420格式，计算公式如下（利用矩阵进行计算）
# Y | 0.299 0.587 0.114 | R 0
# U = | -0.169 -0.331 0.5 | * G + 128
# V | 0.5 -0.419 -0.081 | B 128
# 输入：图片的RGB矩阵
# 输出：二次采样后的YUV格式的Y、U、V
def rgb2yuv420(Image):
# 创建Y、U、V三个数组用于存储RGB->YUV420变换后的值
y = ndarray(shape=(shape(Image)[0], shape(Image)[1]), dtype='float')
u = ndarray(shape=(shape(Image)[0], shape(Image)[1]), dtype='float')
v = ndarray(shape=(shape(Image)[0], shape(Image)[1]), dtype='float')
# 创建跟上面RGB->YUV变换所需要的变换矩阵
matrix_ = matrix([[0.299, 0.587, 0.114], [-0.169, -0.331, 0.5], [0.5, -0.419, -0.081]])
for i in range(shape(Image)[0]):
for j in range(shape(Image)[1]):
image_matrix = matrix(Image[i,j])
# 对每一个像素进行变换
y[i, j], u[i, j], v[i, j] = matrix_*transpose(image_matrix) + matrix([[0],[128],[128]])
# 声明一个等于原来四分之一的数组用与存储二次采样后的U、V的元素
qut_u = ndarray(shape=(shape(Image)[0]//2, shape(Image)[1]//2), dtype='float')
qut_v = ndarray(shape=(shape(Image)[0]//2, shape(Image)[1]//2), dtype='float')
for i in range(shape(qut_u)[0]):
for j in range(shape(qut_v)[1]):
# 二次采样，也就是2✖2的像素只取一个像素点的U、V值，这里统一取第一个点的值
qut_u[i, j] = u[i * 2, j * 2]
qut_v[i, j] = v[i * 2, j * 2]
return y, qut_u, qut_v- DCT变换：
DCT是一种广泛应用的变换编码方法，能够以数据无关的方式解除输入信号之间的相关性。可以将时域转化为频域从而使得原本复杂无关的的数据变得简单。
DCT的变换公式如下：
正变换：
当M与N相等的时候，就可以转变为矩阵运算，如下：
F’= AFAT
代码如下：
# 对矩阵进行分割同时进行dct变换以节省时间
# 输入：需要分割的矩阵与分割成的块的大小，默认为8*8
# 输出：按从左到右从上到下的顺序存储的一系列8*8的已经进行dct变换的块的列表
def split_matrix_and_dct(w_array, block_size=8):
# 声明空的列表
w_arrays = []
width_ = len(w_array[0])
height_ = len(w_array)
# 分成的块的数目等于 （height_/block_size）*（width_/block_size）
for i in range(int(height_/block_size)):
for j in range(int(width_/block_size)):
# 用临时变量存储提取后的8*8的块
temp = w_array[i * block_size: i * block_size + block_size, j * block_size :j * block_size + block_size]
# 在提取的时候顺便进行dct变换，提高效率
m_temp = dct(temp)
# 将dct变换后的数据存储列表中
w_arrays.append(m_temp)
return w_arrays
# DCT变化 利用矩阵进行dct2变换 T=A*F*AT（AT是A的转置）
# 输入：需要进行dct变换的矩阵
# 输出：dct变换后的矩阵
def dct(matrix_):
# 获取矩阵的大小（认为长与宽是一样的）
height_, width_ = matrix_.shape
A = zeros([height_, width_])
N = width_
A[0, :] = sqrt(1/N)
for i in range(1, height_):
for j in range(width_):
A[i, j] = cos(pi * i * (2 * j + 1)/(2 * N))*sqrt(2/N)
# 计算A*F
dct1 = dot(A, matrix_)
# 计算A*F*AT
dct_result = dot(dct1, transpose(A))
# 返回dct变换后的结果
return dct_result- 量化：
在前面的DCT变换后的数据的基础上，将矩阵除以相应的量化矩阵就可以了。（注意有色度与亮度两种量化矩阵），代码如下：
量化矩阵：
Qy = array([[16, 11, 10, 16, 24, 40, 51, 61],
[12, 12, 14, 19, 26, 58, 60, 55],
[14, 13, 16, 24, 40, 57, 69, 56],
[14, 17, 22, 29, 51, 87, 80, 62],
[18, 22, 37, 56, 68, 109, 103, 77],
[24, 35, 55, 64, 81, 104, 113, 92],
[49, 64, 78, 87, 103, 121, 120, 101],
[72, 92, 95, 98, 112, 100, 103, 99]])
Qi = array([[17, 18, 24, 47, 99, 99, 99, 99],
[18, 21, 26, 66, 99, 99, 99, 99],
[24, 26, 56, 99, 99, 99, 99, 99],
[47, 66, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99]])- 量化函数：
# 对dct变换后的矩阵进行量化
# 输入：dct变换后的矩阵的列表 以及 tag标签，tag为‘Y'表示亮度Y的量化，其他则是色度的量化
# 输出：量化的dct变换后的矩阵的列表
def quantization(arr, tag):
# 获取待处理的列表的长度
s = len(arr)
# 判断标签，若为’Y‘则表示为亮度的量化
if tag == "Y":
# 声明空列表
result = []
for i in range(s):
# 除以亮度量表并取整
temp = (arr[i]/Qy).astype("int")
# 加入结果中
result.append(temp)
return result
# 与上面的类似，只是除法的时候采用的是色度量化表
else:
result = []
for i in range(s):
temp = (arr[i] / Qi).astype("int")
result.append(temp)
return result
- AC的RLC编码：
AC，是指不包含量化后的矩阵的M[0][0]的其余数值。
AC的ZigZag编码：也就是Z字形编码，主要是需要发现规律，也就是在碰壁的时候选择有四种情况：
-
- 当是第一行时，则向前进一格后选择往左下方移动
- 当是第一列时，则向下进一格后选择往右上方移动
- 当是最后一行时，则向前进一格后选择往后上方移动
- 当时最后一列时，则向下进一格后选择往左下方移动
代码实现如下：
# AC
# Zig-Zag编码
# 输入：一个矩阵（不是列表）
# 输出：Z编码
def zig_zag(arr):
row, col = arr.shape
num = []
# 声明每一条对角线上的元素的数目
length = [1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1]
# 声明每一个边界条件时数组元素的数目（可以等同于上面length的前i项累加）
l_help = [1, 3, 6, 10, 15, 21, 28, 36, 43, 49, 54, 58, 61, 63, 64]
# 记录判断到length与l_length的位置
l_index = 0
num.append(arr[0][0])
# 用于判断是往右上方还是左下方
up = 0
down = 1
# 用于记录当前遍历的位置的横纵坐标
r = 0
c = 0
# 因为到一半之后规则变换，所以先处理前面的，后面的类似
while len(num) <= 35:
# num的长度与l_help对应的index的长度相等时候就是到达了边界，需要进行判断进行方向的变换
if len(num) == l_help[l_index]:
# 下标前进1
l_index = l_index + 1
# 若是第一行，则说明应该往左下方前进
if r == 0:
down = 1
up = 0
c = c + 1
# 其余情况也就是第一列，说明应该往右上方前进
else:
up = 1
down = 0
r = r + 1
if down == 1:
# 将length中对应的个数的元素添加进num数组中
for i in range(length[l_index] + 1):
if c - i < 0:
r = c
c = 0
break
num.append(arr[r + i][c - i])
elif up == 1:
for i in range(length[l_index] + 1):
if r - i < 0:
c = r
r = 0
break
num.append(arr[r - i][c + i])
# 与上面的类似，只是规则发生了一点点变换
while len(num) <= 63:
if len(num) == l_help[l_index]:
l_index = l_index + 1
if c == col - 1:
down = 1
up = 0
r = r + 1
else:
up = 1
down = 0
c = c + 1
if down == 1:
for i in range(length[l_index] + 1):
if c - r - i < 0:
c, r = r, c
break
num.append(arr[r + i][c - i])
elif up == 1:
for i in range(length[l_index] + 1):
if r - c - i < 0:
c, r = r, c
break
num.append(arr[r - i][c + i])
return numAC的RLC：RLC，游程编码，每一个数值由（a，b）的格式组成，其中b是一个代表AC中出现的非零的一个元素，而a表示的是自上一个非零元素以来的0的数目。统计完成之后以（0，0）结尾。
实现的代码如下：
# 游长编码
# 输入：一个数组
# 输出：一个由(a,b)组成的列表，其中a是b前面据上一个非零元素之间0的数目，b是一个非零元素
def rlc(arr):
# 空列表
pairs = []
# 统计非零元素的数目
count = 0
for i in range(len(arr)):
# 当不等于0的时候将元素与统计的非零元素的数目加入列表中
if arr[i] != 0:
pairs.append([count, arr[i]])
# 重新置为0
count = 0
else:
# 当等于0的时候count++
count = count + 1
# 加入最后一个元素，也就是（0，0）
pairs.append([0, 0])
# 返回列表
return pairs- DC的DPCM编码：
DC，是指所有的8×8的矩阵的M[0][0]所构成的新的数组。
DPCM，也就是差分脉冲调制编码。因为DC所构成的数组中相邻的两个点之间的差别不会很大，因为DC代表的是前面的8×8矩阵中原值的平均值。所以计算的规则如下，创建一个数组arr，其中第一个元素也就是第一个dc值，第i个元素则是dc[i] - dc[i – 1]，也就是相邻的两个dc值之间的差值。
实现的代码如下：
# 获取DC编码
# 输入：一个列表，因为dc是每个方块的第一个元素
# 输出：一个包含所有dc值的数组
def getdc(arr):
# 获取列表的长度并创建一个矩阵
dcs = zeros([len(arr)])
# 将每一个矩阵的第一个元素加入dcs中
for i in range(len(arr)):
dcs[i] = arr[i][0][0]
return dcs
# DPCM编码
# 输入：dc数组
# 输出：根据dpcm进行编码后的数据
def dpcm(arr):
dcs = getdc(arr)
result = zeros(len(arr))
# 结果数组的第一个元素是dc数组的第一个元素
result[0] = dcs[0]
for i in range(len(dcs) - 1):
# 结果数组的第i个元素等于dc数组的第i个元素 - dc数组的第i - 1个元素
result[i + 1] = dcs[i + 1] - dcs[i]
return result- 熵编码前的准备步骤：
DC的VLI（中间格式）编码：为了进一步节省空间，JPEG并不直接保存数据的具体数值，而是将数据按照位数分为16组，保存在表里面。这样就是所谓的变长整数编码VLI。对应的有VLI编码表。转换的规则如下：根据前面的DPCM得到后的数据，比如5，通过VLI表可以看到属于第三组，则表示为（3，5）
代码如下：
# DC VIL中间格式
# 输入：dc的dpcm编码
# 输出：一个[:][2]的矩阵，根据VLI表格将原本[a]扩展为[b, a]，其中b是根据VLI得到的分组（同时也就是bit长度）
def dc_vil(arr):
# 简化VLI表格，下标表示分组
dc_help = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768]
# 声明中间格式矩阵
dcvil = zeros([len(arr), 2])
# 分组
dc_index = 0
for i in range(len(arr)):
# 寻找分组，也就是第一个小于给定的a的dc_help的下标
while abs(arr[i]) >= dc_help[dc_index]:
dc_index = dc_index + 1
# 存储中间格式
dcvil[i] = [dc_index, arr[i]]
dc_index = 0
# 返回结果
return dcvil.astype("int")AC的VLI（中间格式）编码：对前面RLC得到的二元组右边的数据就行如上面的变换变成三元组就可以了。比如（1，5）变为（1，3，5）
代码如下：
# AC VIL中间格式
# 输入：ac的rlc编码
# 输出：一个[:][3]的矩阵，根据VLI表格将原本[a, b]扩展为[a, c, b]，其中b是b根据VLI得到的分组（同时也就是bit长度）
def ac_vil(arr):
# 简化VLI表格，下标表示分组
ac_help = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768]
# 声明中间格式矩阵
acvil = zeros([len(arr), 3])
# 分组
ac_index = 0
for i in range(len(arr)):
# 寻找分组，也就是第一个小于给定的b的ac_help的下标
while abs(arr[i][1]) >= ac_help[ac_index]:
ac_index = ac_index + 1
# 存储中间格式
acvil[i] = [arr[i][0], ac_index, arr[i][1]]
ac_index = 0
# 返回结果
return acvil.astype("int")-
- 熵编码：在得到DC与AC的中间格式之后为了进一步的图像压缩，可以采用哈夫曼编码进行压缩。其中哈夫曼编码时DC与AC系数分别采用不同的编码表，对于色度与亮度也是这样的。因此需要四张编码表才能够完成熵编码，在JPEG标准的有缺省的编码表，在此为了方便所以直接使用所给的编码表。
哈夫曼编码表见：https://www.cnblogs.com/buaaxhzh/p/9119870.html
AC的哈夫曼编码：
上面的例子中，根据（1，3）查表可以得到对应的哈夫曼编码值，而5则用VLI表进行编码，所以这一个元组的长度也就是两者的和。
代码实现如下：
# AC 哈夫曼编码获取 编码的大小
# 输入：ac的中间格式 以及 tag标签，tag为‘Y'表示亮度Y，其他则是色度
# 输出：哈夫曼编码后的长度
def ac_huffman(arr, tag):
# 如果标签为Y则是使用ac的亮度量化表
if tag == "Y":
l = 0
for i in range(len(arr)):
for j in range(len(ACy)):
# 寻找亮度量化表中对应的数值，转化为str判断是否相等
if ACy[j][0] == str(arr[i][0]) and ACi[j][1] == str(arr[i][1]):
# 将VLI的长度与huffman的长度一起存储起来
l = l + int(ACy[j][2]) + arr[i][1]
break
return l
else:
l = 0
for i in range(len(arr)):
for j in range(len(ACi)):
# 寻找色度量化表中对应的数值，转化为str判断是否相等
if ACi[j][0] == str(arr[i][0]) and ACi[j][1] == str(arr[i][1]):
# 将VLI的长度与huffman的长度一起存储起来
l = l + int(ACi[j][2]) + arr[i][1]
break
return lDC的哈夫曼编码：
上面的例子中，用3进行查表得到对应的哈夫曼编码值，而5的利用VLI表进行编码，同上。
代码实现如下：
# DC 哈夫曼编码获取 编码的大小
# 输入：dc的中间格式 以及 tag标签，tag为‘Y'表示亮度Y，其他则是色度
# 输出：哈夫曼编码后的长度
def dc_huffman(arr, tag):
# 如果标签为Y则是使用dc的亮度量化表
if tag == "Y":
# 存储长度
l = 0
for i in range(len(arr)):
for j in range(len(DCy)):
# 寻找亮度量化表中对应的数值，转化为str判断是否相等
if DCy[j][0] == str(arr[i][0]):
# 将VLI的长度与huffman的长度一起存储起来
l = l + int(DCy[j][1]) + arr[i][0]
break
return l
else:
# 存储长度
l = 0
for i in range(len(arr)):
for j in range(len(DCi)):
# 寻找色度量化表中对应的数值，转化为str判断是否相等
if DCi[j][0] == str(arr[i][0]):
# 将VLI的长度与huffman的长度一起存储起来
l = l + int(DCi[j][1]) + arr[i][0]
break
return l- 计算压缩后的数据所需要的比特数目
需要AC的列表计算函数，因为一张图片只有一个dc，所以dc所要的bit长度在上面就可以计算出来了，但是一张图片有多个ac块，所以有：
# 计算所有的AC编码大小
# 输入ac块的列表 与tag标签，tag为‘Y'表示亮度Y，其他则是色度
# 输出：经过z编码、rlc编码、中间格式、与haffman编码之后的总长度
def ac_generator(arr, tag):
if tag == "Y":
l = 0
for i in range(len(arr)):
# 函数嵌套
l = l + ac_huffman(ac_vil(rlc(zig_zag(arr[i]))), "Y")
return l
else:
l = 0
for i in range(len(arr)):
# 函数嵌套
l = l + ac_huffman(ac_vil(rlc(zig_zag(arr[i]))), "I")
return l以下为JPEG的整个压缩过程：
# 主函数
def main():
global width
global height
im = Image.open("动物卡通图片.jpg")
matrix = array(im)
width = len(matrix[0])
height = len(matrix)
# 进行图片的扩展
im = spread_image(im)
matrix = array(im)
width = len(matrix[0])
height = len(matrix)
# 进行图片的颜色空间变换以及二次采样
y, u, v = rgb2yuv420(array(im))
# 进行各个通道上的分块
y_dct = split_matrix_and_dct(y, 8)
u_dct = split_matrix_and_dct(u, 8)
v_dct = split_matrix_and_dct(v, 8)
# 进行各个通道上的量化
y_q = quantization(y_dct, "Y")
u_q = quantization(u_dct, "I")
v_q = quantization(v_dct, "I")
# 计算总的dc编码所需要的bit位数
dc_length = dc_huffman(dc_vil(dpcm(y_q)), "Y") + dc_huffman(dc_vil(dpcm(u_q)), "I") + dc_huffman(dc_vil(dpcm(v_q)), "I")
print(dc_length)
# 计算总的ac编码所需要的bit位数
ac_length = ac_generator(y_q, "Y") + ac_generator(u_q, "I") + ac_generator(v_q, "I")
print(ac_length)显示的结果如下：（第一个为dc的编码所需的bit数，第二个为ac的编码所需的bit数）
解码过程：
- 逆量化：逆量化的过程实际比较简单，也就是在前面量化后的数据（因为保存为整形的，所以会有损失）重新点乘量化矩阵后即可得到逆量化的数据。
# 逆量化
def inv_quantization(arr, tag):
# 获取待处理的列表的长度
s = len(arr)
# 判断标签，若为’Y‘则表示为亮度的逆量化
if tag == "Y":
# 声明空列表
result = []
for i in range(s):
# 乘以亮度量表并取整
temp = (arr[i] * Qy).astype("int")
# 加入结果中
result.append(temp)
return result
# 与上面的类似，只是乘法的时候采用的是色度量化表
else:
result = []
for i in range(s):
temp = (arr[i] * Qi).astype("int")
result.append(temp)
return result-
- 逆DCT（也就是IDCT变换）：
IDCT变换与前面的DCT变换一样，可以使用矩阵进行计算，因为A是一个正交矩阵，所以因为逆变换也就可以表示为。
代码如下：
# IDCT变化 与上面的dct变换类似。也就是将dct变换后的矩阵f，进行T = AT * f * A的变换，其中A与上面的A相同
# 输入：需要进行idct变换的矩阵
# 输出：idct变换后的矩阵
def idct(matrix_):
# 获取矩阵的大小（认为长与宽是一样的）
height_, width_ = matrix_.shape
# 声明A矩阵
A = zeros([height_, width_])
N = width_
A[0, :] = sqrt(1 / N)
# 构造A矩阵
for i in range(1, height_):
for j in range(width_):
A[i, j] = cos(pi * i * (2 * j + 1) / (2 * N)) * sqrt(2 / N)
# 计算AT*f
idct1 = dot(transpose(A), matrix_)
# 计算AT*f*A
idct_result = dot(idct1, A)
# 返回idct变换后的结果
return idct_result- YUV通道扩展合并：
将8×8的数据列表整合成一个完整的矩阵。全局变量height与width记录了图片扩展之后的大小。其中Y矩阵的大小与图片一样，U、V矩阵因为进行了二次采样，因此大小是图片的四分之一。然后根据自左向右，自上向下的规则进行合并即可。
# 将idct后三通道的数据进行整合
# 输入：需要合并的通道arr， 标签tag（Y代表y通道，其余字母表示u，v通道），block_size获取arr中每一个元素的大小
# 输出：合并后的通道的矩阵mat1
def generate_m(arr, tag, block_size=8):
if tag == "Y":
# height、width是全局变量，也就是原来图片的长高
mat1 = zeros([height, width])
for i in range(height//block_size):
for j in range(width//block_size):
print(i * width // block_size + j)
# 在整合之前进行idct变换以缩短时间提高效率
mat1[i * block_size: i * block_size + block_size, j * block_size:j * block_size + block_size] = idct(
arr[i * width // block_size + j])
return mat1
else:
# 因为U、V通道进行了二次采样，是原来大小的1/4
mat1 = zeros([height//2, width//2])
for i in range(height // block_size//2):
for j in range(width // block_size//2):
# 在整合之前进行idct变换以缩短时间提高效率
mat1[i * block_size: i * block_size + block_size, j * block_size:j * block_size + block_size] = idct(
arr[i * width // block_size // 2 + j])
return mat1- YUV到RGB的颜色空间转化：
这跟前面的一系列逆变换是一样的，也是利用公式的逆推导进行计算即可
# 将YUV转化为RGB格式，运用上面的公式进行逆变换，上面的矩阵取逆
# R | 0.299 0.587 0.114 | ' (Y 0)
# G = | -0.169 -0.331 0.5 | * (U - 128)
# B | 0.5 -0.419 -0.081 | (V 128)
# 输入：图片的RGB矩阵
# 输出：二次采样后的YUV格式的Y、U、V
def yuv2rgb(mat1):
# 创建跟上面RGB->YUV变换所需要的变换矩阵
print(mat1.shape)
matrix_ = matrix([[0.299, 0.587, 0.114], [-0.169, -0.331, 0.5], [0.5, -0.419, -0.081]])
# 创建RGB空间矩阵
rgb_image = zeros([len(mat1), len(mat1[0]), 3])
for i in range(len(mat1)):
for j in range(len(mat1[0])):
# 上面的公式表示
rgb_image[i][j] = squeeze((mat(matrix_).I * (matrix([[mat1[i][j][0]], [mat1[i][j][1]], [mat1[i][j][2]]]) - matrix([[0], [128], [128]]))).transpose())
# 检测表换后的数值是否在0 - 255之间，若过大或过小，则置为最大或最小值
if rgb_image[i][j][0] > 255:
rgb_image[i][j][0] = 255
if rgb_image[i][j][1] > 255:
rgb_image[i][j][1] = 255
if rgb_image[i][j][2] > 255:
rgb_image[i][j][2] = 255
if rgb_image[i][j][0] < 0:
rgb_image[i][j][0] = 0
if rgb_image[i][j][1] < 0:
rgb_image[i][j][1] = 0
if rgb_image[i][j][2] < 0:
rgb_image[i][j][2] = 0
return rgb_image- 完整的解码过程代码示例：
# 三通道的逆量化
i_y_q = inv_quantization(y_q, "Y")
i_u_q = inv_quantization(u_q, "I")
i_v_q = inv_quantization(v_q, "I")
# 三通道的逆dct变换与整合
i_y = generate_m(i_y_q, "Y")
i_u = generate_m(i_u_q, "I")
i_v = generate_m(i_v_q, "I")
# 利用三通道生成图片并进行颜色空间的转化
m0 = (yuv2rgb(generate_matrix(i_y, i_u, i_v))).astype(int)
# 修改为image的格式
m0 = uint8(m0)
# 选取图片原来的大小的矩阵（因为有一部分是扩展的）
result = m0[0:height1, 0:width1, :]
# 利用矩阵生成图片
im2 = Image.fromarray(result)
# 存为bmp格式以便后面计算误差
im2.save("JPEG.bmp")结果对比分析：
视觉效果的对比：
- 采用局部对比（100*100）
-
图片
原图（局部）
JPEG
GIF
动物卡通图片
动物照片
- 视觉参数计算：（代码在后面一点给出）
MSE（图片均方差）：是原图片与生成的图片差值矩阵每一个像素点RGB三通道的平方和，数值越大，说明效果越差。
| 图片 | JPEG | GIF |
| 动物卡通图片 | 0.015710 | 1.676224 |
| 动物照片 | 0.364148 | 1.179221 |
PSNR（峰值信噪比）：等于10log10(MAX2MSE)，表示的是信号与噪声的比值，比值越大，说明效果越好。若为45-55db，则通常认为良好，若为60，则可以认为是优良。
| 图片 | JPEG | GIF |
| 61.397952 | 41.116270 | |
| 动物照片 | 47.746817 | 42.643638 |
SSIM（结构相似性）：是一种衡量两幅图片相似度的指标，计算表达式如下：
SSIMx,y= (2μxμy+c1)(2σxy+c2)(μx2+μy2+c1)(σx2+σy2+c2 )
其中μ是平均值，σ是标准差，c1=(k1L)2, c2=(k2L)2 是用来维持稳定的常数。L是像素值的动态范围。k1=0.01, k2=0.03。结构相似性的范围是0到1，当两张图一模一样的时候SSIM的值为1.
| 图片 | JPEG | GIF |
| 动物卡通图片 | [0.98862801 0.99188251 0.99456432] | [0.95330913 0.9617353 0.98109455] |
| 动物照片 | [0.96031301 0.95045469 0.95294486] | [0.97184365 0.96874164 0.96583153] |
所以通过以上的MSE、PSNR与SSIM三个参数可以看出JPEG压缩的视觉效果基本碾压GIF。
-
- 代码实现如下：
from numpy import *
from PIL import Image
from math import *
def main():
# 打开BMP24位图片
bmp_image = Image.open("./Image/BMP.bmp")
bmp_matrix = array(bmp_image)
# 打开压缩后的JPEG图片
jpeg_image = Image.open("./Image/JPEG.jpg")
jpeg_matrix = array(jpeg_image)
# 打开压缩后的GIF图片
gif_image = Image.open("./Image/GIF.gif")
# 利用getpalette获取256色表填充GIF图片使其变成[:][:][3]维的图片以便进行计算
gif_matrix = gif2rgb(array(gif_image), gif_image.getpalette())
# 获取BMP与JPEG图片像素值的差值
sub_ = bmp_matrix - jpeg_matrix
# 计算矩阵的f范数，也就是每个点的平方之和
sub_jpeg_fro = calculate(sub_)
# 计算MSE均方差 等于f范数除以矩阵中元素的个数
mse_jpeg = sub_jpeg_fro/(len(jpeg_matrix[0])*len(jpeg_matrix)*3)
# 输出JPEG的MSE 均方差
print("the MSE of the JPEG is %f" % mse_jpeg)
# 判断差值矩阵的F范数是否为0，若是则说明无损，也就是PSNR为inf
if sub_jpeg_fro == 0:
print("the SNR of the JPEG is inf")
# 否则进行计算
else:
# 计算PSNR 峰值信噪比 10*log10(255*255*矩阵的大小*每个像素有RGB三个/F范数)
snr_jpeg = 10 * log10(255**2*len(jpeg_matrix[0])*len(jpeg_matrix)/sub_jpeg_fro)
# 输出JPEG的PSNR峰值信噪比
print("the SNR of the JPEG is %f" % snr_jpeg)
# 调用ssim函数输出SSIM的结果
print("the SSIM of the JPEG is " ,cal_ssim(bmp_matrix, jpeg_matrix))
# 与上面JPEG相似，此处不做累述
sub_gif_fro = calculate((bmp_matrix - gif_matrix))
mse_gif = sub_gif_fro/(len(gif_matrix[0])*len(gif_matrix)*3)
print("the MSE of the GIF is %f" % mse_gif)
if sub_gif_fro == 0:
print("the SNR of the GIF is inf")
else:
snr_gif = 10 * log10(255**2*len(gif_matrix[0])*len(gif_matrix)/sub_gif_fro)
print("the SNR of the GIF is %f" % snr_gif)
print("the SSIM of the GIF is" ,cal_ssim(bmp_matrix, gif_matrix))
# 用于计算矩阵的F范数
# 输入：矩阵
# 输出：F范数
def calculate(mat_3d):
# 提取各个通道的数值
mat0 = mat_3d[:][:][0]
mat1 = mat_3d[:][:][1]
mat2 = mat_3d[:][:][2]
# 利用 mat*mat得到一个每个位置上元素的平方的矩阵，再利用sum进行求和
sum_ = sum(mat0*mat0) + sum(mat1*mat1) + sum(mat2*mat2)
return sum_
# 利用颜色查找表来扩充GIF图片来使其变成[:][:][3]的矩阵以便后面计算
# 输入：GIF原矩阵 颜色查找表
# 输出：[:][:][3]的矩阵
def gif2rgb(mat_d, row_palette):
# 利用get_palette对获取的原颜色查找表进行处理
palette = get_palette(row_palette)
# 声明一个[:][:][3]的矩阵
mat3_d = zeros([len(mat_d), len(mat_d[0]), 3])
for i in range(len(mat_d)):
for j in range(len(mat_d[0])):
# 因为GIF每个像素点存的是颜色查找表的索引，因此利用索引找到对应的颜色即可
mat3_d[i][j] = array(palette[mat_d[i][j]])
return mat3_d
# 原颜色查找表的格式是[;][1]的列表，需要将其转化为[:][3]的矩阵
# 输入：颜色查找表列表
# 输出：颜色查找表矩阵
def get_palette(arr):
s = len(arr)
# 声明一个[:][3]的矩阵
palettes = zeros([s//3, 3])
for i in range(s//3):
# 每一行的赋值
palettes[i] = array([arr[3*i], arr[3 * i + 1], arr[3 * i + 2]])
return palettes
# 计算SSIM参数
# 输入：两个矩阵，以及SSIM所需要的常熟 k1,k2,l,通常为0.01，0.03，255
# 输出：一个三个元素的数组，对应每一个通道上的SSIM值
def cal_ssim(mat1, mat2, k1 = 0.01, k2 = 0.03, l = 255):
# 计算每个矩阵的每个通道上平均值ux, uy
mean1 = array([mat1[:][:][0].mean(), mat1[:][:][1].mean(), mat1[:][:][2].mean()])
mean2 = array([mat2[:][:][0].mean(), mat2[:][:][1].mean(), mat2[:][:][2].mean()])
# 调用函数计算每个矩阵的方差（因为numpy中的var函数比较奇怪，所以自己实现了一个）
# sigmax^2, sigmay^2
var1 = cal_var(mat1, mean1)
var2 = cal_var(mat2, mean2)
# 用于存储x、y协方差的每一个元素值，也就是(X(i, j) - ux)*(Y(i, j) - uy)
sigma_ = []
for i in range(len(mat1)):
for j in range(len(mat1[0])):
# 计算并存储(X(i, j) - ux)*(Y(i, j) - uy)
sigma_.append((mat1[i][j] - mean1) * (mat2[i][j] - mean2))
# 计算和用于计算sigmaxy
sigma_sum = sum(sigma_, axis=0)
sigma = sigma_sum/(len(mat1) * len(mat1[0]) - 1)
# 规定的常数计算
c1 = (k1 * l)**2
c2 = (k2 * l)**2
# 扩展为三个元素的数组以便后面的运算
c1_arr = array([c1, c1, c1])
c2_arr = array([c2, c2, c2])
# (2ux*uy + c1)*(2sigmaxy + c2)
# 计算SSIM参数， SSIM(X, Y) = --------------------------------------------
# (ux^2 + uy^2 + c1)(sigmax^2 + sigmay^2 + c2)
ssim = (2 * mean1 * mean2 + c1_arr) * (2 * sigma + c2_arr) /((mean1 * mean1 + mean2 * mean2 + c1_arr) * (var1 + var2 + c2_arr))
return ssim
# 计算方差
# 输入：矩阵、矩阵的平均值
# 输出：矩阵的方差(每个通道上)
def cal_var(mat, mean):
sigma_sum = [0, 0, 0]
for i in range(len(mat)):
for j in range(len(mat[0])):
sigma_sum += (mat[i][j] - mean)**2
res = sigma_sum/(len(mat)*len(mat[0]) - 1)
return array(res)
if __name__ == '__main__':
main()- 压缩效果的对比：
首先将图片另存为BMP格式，与自己生成之后的JPEG与GIF的大小进行对比并计算。压缩率越小，说明压缩的程度越高。
| 图片 | BMP | JPEG 熵编码长度 | 熵编码后的大小/bmp的大小 | GIF | GIF压缩率 |
| 动物卡通图片 | 2088.96KB | 95.97KB | 4.58% | 347KB | 16.61% |
| 动物照片 | 2037.76KB | 115.8KB | 5.66% | 410KB | 20.12% |
上面的结果显示JPEG的压缩程度会更高。
三、实验心得
- 平时上课没怎么听懂的知识通过这一次作业搞懂了很多，因为JPEG压缩涉及的东西很多，比如DCT离散余弦变换，还有DPCM编码、RLC编码、颜色空间转化与二次采样和哈夫曼编码。其中二次采样的步骤比较坑，因为YUV420当时听课的时候就觉得比较玄学，而且在网上查找资料的时候一开始也没有找到，找了很久之后才知道YUV420与YUV411的区别只是一小点点，但是却有着不同优势与劣势。其余的部分相对比较简单，只要照着公式实现就好了。
- 这一次实验自己写了1000多行代码（不包含常量：量化表与哈夫曼表），可以说是印象十分深刻了。
670
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
