使用 PCA 实现图像压缩

使用 PCA 实现图像压缩

我们知道,一张灰度照片在电脑里就是一张数字表格,H 行 W 列,每个格子存一个 0 到 255 的亮度值,彩色图再叠上红、绿、蓝三张这样的表格。可仔细想,表格里大部分数字是重复的,天空那一片几乎一个颜色,墙面那一块明暗缓慢变化,相邻格子的数值高度相关。真正在画面里制造变化的信息量,远没有像素总数那么多。主成分分析 PCA 就是冲着这件事来的,它找出数据里变化最大的几个方向,用每个点在这几个方向上的新坐标,代替原来成百上千个坐标,数据就压下来了。

这一节我们先从方差和特征值讲清 PCA 到底在找什么方向,再手算一个 5 行 2 列的小矩阵,把中心化、协方差矩阵、特征值、投影、重建一步步走完,最后回到图像压缩的完整代码。读之前需要知道矩阵乘法怎么算,知道特征值和特征向量大致是什么意思,用到时我会先把概念说成人话。

1. 图像即矩阵

我们先把图像翻译成 PCA 能吃的形式。一张灰度图是 H 行 W 列的矩阵 A。做 PCA 时,约定一行是一个样本,一列是一个特征,也就是把这张图看成 H 个样本,每个样本是一个 W 维向量,这个向量记录了某一行像素从左到右的亮度。彩色图有三个通道,就把红、绿、蓝三张表分别当独立的矩阵处理。

为什么这样压得动?因为相邻两行的亮度曲线往往长得差不多,这一整行 W 个数,其实是被少数几种基本曲线拼出来的。PCA 要做的,就是自动找出这几种基本曲线,然后每一行的亮度只要用它在这几条曲线上的投影系数来表示就够了。原来存 H 乘 W 个数,现在存 H 乘 k 个系数,k 远小于 W,体积就降下来了。

2. 方差定方向

那几条基本曲线该怎么找?PCA 给的判据是保留最多的方差。我们把每个样本想象成二维平面上的一个点,压缩到一维就是把这些点投影到某条直线上,只保留每个点在这条线上的坐标。问题是这条线朝哪个方向摆。摆得好,投影后点还能铺开,原来谁是谁分得清。摆得不好,所有点投影挤成一堆,原来的差别全丢了。

所以选方向的标准很明确,让投影之后的方差最大。方差就是这些投影点散开的程度,散得越开,说明这个方向上保留的变化越多,信息丢得越少。怎么自动找这条方差最大的方向?先把数据中心化,每个特征减去自己的均值,让点云居中在原点。再算协方差矩阵,它是一个方阵,对角线上是每个特征自己的方差,非对角线上是两个特征一起变动的程度。这个矩阵的特征向量,正好指向数据云散开的那些主轴方向,对应的特征值,就是那个方向上的方差大小。

一句话记住,特征值越大的方向越重要,它扛着越多的方差。把特征值从大到小排,前面几个对应的特征向量,就是要保留的主成分。后面小特征值对应的方向,方差小,变化微不足道,丢掉也不太影响画面。

3. 小矩阵算例

我们拿一个最小的例子把流程走一遍。假设手头有 5 个样本,每个样本 2 个特征,排成 5 行 2 列的矩阵 X,想把它从 2 维压到 1 维。原始数据长这样:

样本特征 1特征 2
122
243
335
465
512

第一步中心化。两个特征各自的均值是 3.2 和 3.4,每个数减去对应均值,得到居中后的点云。第二步算协方差矩阵。两个特征之间算出来是这样,主对角分别是 3.7 和 2.3,非对角是 2.15,协方差矩阵是个 2 行 2 列的对称方阵。第三步求特征值和特征向量。算出来两个特征值,大的那个约 5.2611,小的那个约 0.7389。大特征值对应的特征向量约是 (0.8092, 0.5875),这就是第一主成分方向。

特征值加起来是 5.2611 加 0.7389,等于 6.0,这就是数据的总方差。第一主成分扛了多少,5.2611 除以 6.0,约 87.68%。第二主成分只扛 0.7389 除以 6.0,约 12.32%。也就是说,把数据压到 1 维,扔掉第二主成分,理论上还能保住八成八的方差。

第四步投影。把每个中心化后的点和这个单位方向向量 (0.8092, 0.5875) 做内积,得到每个点的一维坐标。第五步如果想看看丢了多少,就用这个一维坐标乘回方向向量,再加回均值,重建出一个 2 维点。原始值和重建值放在一起对照:

样本原始 (特征1, 特征2)重建 (特征1, 特征2)
1(2, 2)(1.75, 2.35)
2(4, 3)(3.53, 3.64)
3(3, 5)(3.83, 3.86)
4(6, 5)(5.79, 5.28)
5(1, 2)(1.09, 1.87)

我们看到重建值和原值很接近,但不是完全相等,丢掉的那点差别,就是第二主成分里那 12.32% 的方差。这个小例子在 numpy 和 sklearn 里都能复现,跑一遍就知道上面的数不是编的:

import numpy as np
from sklearn.decomposition import PCA

X = np.array([[2.0, 2.0], [4.0, 3.0], [3.0, 5.0], [6.0, 5.0], [1.0, 2.0]])

p = PCA(n_components=1)
Z = p.fit_transform(X)

print(p.components_)                 # [[0.8092, 0.5875]]
print(p.explained_variance_ratio_)   # [0.8768]
print(np.round(Z, 4).ravel())       # [-1.7936  0.4123  0.7782  3.2058 -2.6028]
print(np.round(p.inverse_transform(Z), 4))

4. 保留多少维

刚才那个例子只有 2 维,一眼就能看到两个特征值。真实图像是几万维,主成分成百上千个,得有个办法自动决定留几个。做法是把特征值从大到小排,逐个累加,看累计方差占比爬到多高。例如想要保住 80% 的信息,就从最大的那个特征值开始往上加,加到累计和超过总方差的 80%,用掉几个特征值就留几个主成分。

sklearn 的 PCA 把这件事做成了参数。传 n_components 等于 0.8,它不是说留 0.8 个分量,而是自动选出最少的 k,让前 k 个主成分的累计方差占比不低于 0.8。我们只说保留多少信息,具体压到几维由数据自己定,图像越有规律,需要的 k 越小。

5. 压缩代码

原理落到图像上,代码其实很短。我们读入图片,像素值除以 255 归一到 0 到 1。彩色图三个通道,对每个通道单独做一次 PCA,保留 80% 方差,fit_transform 把矩阵压成低维,inverse_transform 再映射回原来的行列数。三个通道拼回去,把数值截断到 0 到 1,就能显示和保存。完整代码如下:

import matplotlib.pyplot as plt
import numpy as np
from sklearn.decomposition import PCA
import os

if __name__ == '__main__':

    _, axes = plt.subplots(1, 5, figsize=(20, 8))

    data = plt.imread('data/img.webp') / 255
    # 显示原图
    axes[0].imshow(data)
    axes[0].set_title('原图')

    new_img = []
    for index in range(3):
        # 保留图像信息
        transfer = PCA(n_components=0.8)
        chanel = data[:, :, index]
        chanel = transfer.fit_transform(chanel)
        chanel = transfer.inverse_transform(chanel)
        new_img.append(chanel)

    # 显示压缩图
    new_img = np.array(new_img).transpose([1, 2, 0])
    # 将值截断到 0-1 之间
    new_img = np.clip(new_img, 0, 1)
    axes[1].imshow(new_img)
    axes[1].set_title('压缩图')

    for index in range(3):
        axes[2 + index].imshow(new_img[:, :, index])
        axes[2 + index].set_title('通道' + str(index + 1))
    plt.show()


    plt.imsave('data/new_img.webp', new_img)

    print('原图大小:', os.path.getsize('data/img.webp'))
    print('压缩图大小:', os.path.getsize('data/new_img.webp'))

跑出来的文件大小如下:

原图大小: 37848
压缩图大小: 18120

我们看到原图在磁盘上占 37848 字节,压缩后存成 new_img.webp 占 18120 字节,18120 除以 37848 约等于 0.48,体积差不多省了一半。要说明一句,这里量的是 webp 文件大小,webp 本身也是有损编码,所以一半的降幅里既有 PCA 丢信息的功劳,也有 webp 编码的功劳,不能全算在 PCA 头上。代码里那句 inverse_transform 别删,PCA 之后的坐标是在新主成分坐标系里的,而图像必须用原来的行列坐标才能显示,这一步就是把数据映射回原空间。压缩图整体看着还行,只是边缘和细纹比原图糊一点,那就是被丢掉的 20% 方差。

6. 常见误区

第一个坑,是以为 PCA 压缩是无损的。它保住的是方差最大的那部分,细节纹理本来方差就小,恰恰被丢了,所以重建图一定比原图糊,这是有损压缩。

第二个坑,是以为 n_components 等于 0.8 就意味着体积压到原来的 80%。它说的是保留 80% 的方差,实际压到几维、体积小多少,要看图像本身的相关性,画面越平省得越多。

第三个坑,是以为彩色图要把 H 乘 W 乘 3 当一个大矩阵整体送进 PCA。上面的代码是三个通道各做各的,红、绿、蓝互不干扰,这是常见也够用的做法。

第四个坑,是搞反样本和特征。sklearn 的 PCA 默认一行是一个样本、一列是一个特征,对图像来说就是 H 行样本、W 维特征,别把行列转置反了,否则主成分的含义就变了。

图像是矩阵,矩阵里相邻行列高度相关,PCA 就靠找方差最大的方向来去掉这种冗余。流程是先中心化,再算协方差矩阵,求特征值特征向量,按特征值大小挑前 k 个方向,把每个样本投影到这 k 个方向上,用新坐标代替原来的高维坐标。k 由累计方差占比自动定,sklearn 里一句 n_components 等于 0.8 就行。压完用 inverse_transform 映射回原空间才能看图,整个过程是有损的,丢的主要是细节纹理。