数据科学与工程算法课程项目-PCA图像压缩
摘要
本项目针对三种类型，每种类型包含100个样本的256px*256px的图片数据集使用PCA方式进行图像压缩。本项目分析了传统PCA方法应用于图像压缩时可能存在的缺陷，并为了克服这些缺陷提出了Multi-Block PCA和2D PCA两种改进后的PCA方法，对这三种PCA方法进行了设计、实现和比较，分析不同算法在不同数据集、不同压缩率下的图像质量和算法效率，探究在不同应用场景下最适合该场景的图像压缩PCA算法。
This project is aimed at three types, each type contains 100 samples of 256px*256px image data set using PCA method for image compression. This project analyzes the possible shortcomings of the traditional PCA method when applied to image compression, and in order to overcome these shortcomings, two improved PCA methods, Multi-Block PCA and 2D PCA, are proposed, and these three PCA methods are implemented and tested. And compare, analyze the image quality and algorithm efficiency of different algorithms under different data sets and different compression ratios, and explore the image compression PCA algorithm that is most suitable for the scene in different application scenarios.
1. 项目概述
在多元统计分析中，主成分分析（Principal components analysis，PCA）是一种统计分析、简化数据集的方法。它利用正交变换来对一系列可能相关的变量的观测值进行线性变换，从而投影为一系列线性不相关变量的值，这些不相关变量称为主成分（Principal Components）。具体地，主成分可以看做一个线性方程，其包含一系列线性系数来指示投影方向。PCA对原始数据的正则化或预处理敏感（相对缩放）。
传统的PCA方法将图像转化为维数为m的向量，生成n行m列的数据矩阵X，包含n个样本，每个样本包含m维属性。对该数据矩阵采用PCA方法进行压缩，选取k个主成分，将每个图片压缩成一个长度为k的向量，并通过一个$k\times m$的矩阵重构图像，起到图像压缩的效果。
但在该项目中将说明，传统PCA方法并不适用于所有应用场合，尤其在图像压缩时，由于图像数据的特点（高维度，维度间存在相关性），传统PCA对于这种特点的数据效果不够好。因此该项目选择了一些改进后的PCA方法代替传统PCA，分析和比较了这些方法之间的异同和优劣。
该项目主要提出了Multi-Block PCA和2D PCA两种方法解决该问题。
2. 问题描述
在本项目中，数据集图片来自飞机、农田、海岸三种类型，每种类型由100张256px*256px的彩色图片组成。我们需要对彩色图片进行PCA压缩，以压缩率和PSNR作为压缩标准，尽可能以较低的压缩率获得较高的PSNR，对于彩色图片，分别对RGB三个通道计算PSNR值并取平均作为评价标准。
数据集：
效果：
为此，该项目根据已有的研究成果，对给定的数据集分别尝试了传统PCA、Multi-Block PCA和2D-PCA三种方法进行图像压缩，并分别对结果和性能进行测试和比较。
3. 方法
在介绍Multi-Block PCA以及2D-PCA之前，首先对传统的PCA方法进行介绍，并解释传统PCA方法在该问题下的不足之处。
Traditional PCA
传统PCA方法步骤如下：
- 计算中心化向量$\overline x = \frac{1}{n}\sum\limits_{i=1}^n x_i$，并去中心化$y_i=x_i-\overline x$
- 构造$d\times n$维矩阵$Y=[y_1,y_2,\cdots,y_n]$，计算$d\times d$维协方差矩阵$C=\frac{YY^T}{n-1}$
- 对$C$矩阵进行正交对角化（或对$Y$矩阵进行SVD分解），假定$\lambda_i$对应的特征向量为$u_i$
- 保留前$k$个大于$\alpha$的特征对应的特征向量
- 降维重建向量$\hat{z_i}=[u_1,\cdots,u_k]^Tx_i$
- 重构原矩阵$\hat {x_i} = [u_1,\cdots, u_k]\hat {z_i} + \overline x$
其压缩率$CR$计算公式为$CR=\frac{k\times n_samples+k\times M\times N}{M\times N\times n_samples}$，其中k为选取的主成分个数，n_samples为样本个数，M、N分别为图象长宽像素数，可以反推得到$k=\frac{(CR\times M\times N\times n_samples)}{N_samples+M\times N}$
传统方法在该项目中主要有两个缺陷：
- 图像为高维数据，在PCA过程中需要计算大小为$d\times d$协方差矩阵，协方差矩阵过大，使得计算的代价较大，同时也可能导致一定的误差
- 将二维图像转化为一维向量时忽略了图像的特性，即维度之间的相关性，例如利用图像的局部性，相邻的像素点通常有较强的相关性，而传统的PCA方法在二维转一维的过程中破坏了这些相关性
为了减少计算代价，可以对数据矩阵进行SVD分解代替对协方差矩阵的特征分解，减少计算量。由线性代数知识可知$YY^T$与$Y^TY$有相同的特征值，那么可以先对$Y^TY$这个大小为$n\times n$的矩阵进行特征值分解，求出特征值后再由幂法近似迭代求出特征值对应的特征向量。
实际上这一步方法很多，也有大量针对该问题的算法研究，如sklearn矩阵分解模块提供的基于矩阵变换的SVD分解等等，但这些算法研究不在本项目的讨论范围之内。
为了解决第二个问题，接下来介绍该项目主要使用的两种改进后的PCA方法。
Multi-Block PCA
Multi-Block PCA首先将原图像首先切分成若干的小块（如切成16*16个16px*16px的小块）。利用图像的局部性和维度之间的相关性，同一块内的数据可能存在一定的相关性，而不同块之间的相关性相对会较小。基于这样一种假设，我们可以对每一块局部的小块独立进行压缩，使得PCA过程仅关注分块内的信息，不会受其他块信息的干扰，从而有效降低压缩率。
如图所示，图片 B. Qiu, V. Prinet, E. Perrier and O. Monga, “Multi-block PCA method for image change detection”
即将数据矩阵切分成如下形式：
其中
$bs$代表块大小(blocksize)，$u=(i-1)bs, v = (j-1)bs$
其压缩率$CR$计算公式为$CR=\frac{k\times n_samples+k\times bs\times bs}{bs\times bs\times n_samples}$，其中k为选取的主成分个数，n_samples为样本个数，可以反推得到$k=\frac{(CR\times bs\times bs\times n_samples)}{100+bs\times bs}$
具体步骤描述如下：
- 确定分块大小blocksize
- 提取分块信息转化为一维向量$x_i$，将所有图片对应位置的分块构建为数据矩阵$X$
- 计算$X$的协方差矩阵$C$
- 对C进行特征值分解，保留前$k$大特征值对应的特征向量，构建特征矩阵$V$
- 计算降维向量$Z=V^T X$
- 重构分块$\hat X = VZ + X_{mean}$
- 重复步骤2-6，计算每个分块的重构分块
- 将所有重构分块拼接在一起组成重构图像
2D PCA
相较于传统PCA，2D PCA不用将图片转换为1D向量，极大的减少了数据的维度。基本思路和PCA相同，也是将数据投影到某组基上，然后使得投影之后数据的协方差矩阵成为对角阵。投影后数据的协方差矩阵为：
$S_x=E[(Y-EY)(Y-EY)^T]=E[AX-E(AX)][AX-E(AX)]^T=E[(A-EA)X][(A-EA)X]^T$
利用$tr(AB)=tr(BA)$公式变形后得到$S_x=X^T[E(A-EA)^T(A-EA)]X$
再定义矩阵$G_t=E(A-EA)^T(A-EA)$，即原数据的协方差矩阵
接着对$G_t$矩阵进行特征值分解，取前k大的特征向量进行数据压缩
其压缩率$CR$计算公式为$CR=\frac{k\times n_samples\times M+k\times M}{M\times N\times n_samples}$，其中k为选取的主成分个数，n_samples为样本个数，可以反推得到$k=\frac{(CR\times N\times M\times n_samples)}{n_samples\times M+M}$
具体步骤如下：
- 求出所有图像的平均矩阵$average = \overline X$
- 根据上式计算$G_t$矩阵，即$G_t=\sum\limits_i(X_i-average)^T(X_i-average)$
- 对$G_t$进行特征值分解，保留前$k$大特征值对应的特征向量，构建特征矩阵$V$
- 计算降维向量$Z=(X-average)V$
- 重构矩阵$\hat X = ZV^T+average$
4. 实验结果
测试block_size的影响
为了进一步比较三种PCA算法，需要先为Multi-Block PCA找到合适的block_size。通过设置不同的block_size参数，计算在该参数下Multi-Block在不同压缩率下的表现情况。可以发现三种数据集（飞机、海岸、农田）图像的复杂程度递减，因此选择压缩难度中等的’beach’作为测试数据集。
测试标准：有效性：计算PSNR，度量原图像和重构图像的误差
其中$MSE=\sum\limits_{i}\sum\limits_{j}\sum\limits_{k}(X[i,j,k]-\hat X[i,j,k])^2, MAXI=255$
高效性：测量程序运行时间
| block_size = 4 | |||
|---|---|---|---|
| CR | 0.25 | 0.5 | 0.75 |
| PSNR | 38.65 | 43.60 | 50.64 |
| block_size = 8 | |||
|---|---|---|---|
| CR | 0.25 | 0.5 | 0.75 |
| PSNR | 37.69 | 43.90 | 49.23 |
| block_size=16 | |||
|---|---|---|---|
| CR | 0.25 | 0.5 | 0.75 |
| PSNR | 36.26 | 42.61 | 48.47 |
| block_size=32 | |||
|---|---|---|---|
| CR | 0.25 | 0.5 | 0.75 |
| PSNR | 34.02 | 40.32 | 47.51 |
| block_size | runtime(s) |
|---|---|
| 4 | 7.85 |
| 8 | 6.21 |
| 16 | 60.49 |
| 32 | 466.69 |
综合比较不同block_size下的重构图像质量和程序运行时间，选择block_size=8作为Multi-Block算法的超参数，此时算法的性能表现和运行效率都比较高
比较算法有效性
基于以上三种PCA方法的实现，设计了多组实验测试方法的有效性和高效性，主要从压缩率和PSNR两个指标作为评价的量化标准。
在保证压缩率分别为5%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90%, 95%的情况下选择一张图片直观展示压缩后的图像情况，并对PSNR进行测试。
为了保证压缩率符合要求，首先根据压缩率反推出对应算法所需要保留的主成分维度k，然后设置参数k，使得对应算法在保留k个特征向量的情况下进行压缩测试，保证实际压缩率尽可能接近期望的压缩率。
直观展示
原图像：
Traditional PCA：
Multi-Block PCA：
2D PCA：
原图像：
Traditional PCA：
Multi-Block PCA：
2D PCA：
原图像：
Traditional PCA：
Multi-Block PCA：
2D PCA：
PSNR
绘制结果如下：
由图像可知，在普遍情况下（5%<=CR<=80%）压缩效果Multi-Block PCA > 2D-PCA > traditional PCA，当压缩率趋向于1时，各种算法表现不稳定，可能出现2D-PCA和传统PCA表现更好的情况。但考虑到这些情况实际意义较小（压缩率过大时相当于没有压缩），因此普遍来说可以认为Multi-Block算法在该数据集下表现最佳。
横向比较三张图，可以发现airplane中三种算法表现差异最大，agricultural表现差异最小。这可以由数据集图像本身的复杂程度解释，airplane复杂程度最高，因此对算法的要求更高，糟糕的算法更容易产生糟糕的结果；而agricultural图像最简单，同一图像内和图像间都变化较小，因此即使使用较差的压缩算法，结果仍然可接受。
比较算法高效性
通过测量算法运行时间来作为算法高效性判断依据。
类似于验证算法有效性的方法，绘制不同算法在不同数据集下以不同压缩率运行所需要的时间
从运行时间的角度来说，三种算法的效率2D-PCA > traditional PCA > Multi-Block PCA。
这与三种算法的实现方式有关，传统PCA和Multi-Block PCA都需要将二维数据转为一维再进行PCA降维，而2D-PCA没有这一步，这使得PCA过程中的维度数量有很大的差别，进而导致算法复杂度不同。粗略分析，对于$m$个$n\times n$的图像，单次2D-PCA的复杂度为$O(mn^2)$，传统PCA和Multi-Block PCA的时间复杂度则远大于该复杂度（具体时间复杂度依赖于PCA算法中所采取的SVD算法复杂度，在此不进行分析）。
5. 结论
根据以上测试结果，从有效性和高效性两个角度进行比较，可以发现传统PCA的确在数据集上的表现不如Multi-Block PCA和2D PCA两种改进后的算法。
比较Multi-Block PCA和2D PCA则可以发现，在压缩率相同的情况下Multi-Block PCA生成的压缩图像质量稍好于2D PCA，但2D PCA在运行时间上明显好于Multi-Block PCA。说明Multi-Block PCA和2D PCA分别适用于不同的应用场景，Multi-Block适用于对压缩质量要求较高，但对压缩时间要求不高的场合；2D PCA适用于对压缩时间要求较高的场合，同时压缩质量也能有一定的保证。
在选择压缩率时，可以根据上面绘制的压缩率-PSNR图进行选择。根据PSNR指标，PSNR高于40dB说明图像质量极好（即非常接近原始图像），在30—40dB通常表示图像质量是好的（即失真可以察觉但可以接受），在20—30dB说明图像质量差，PSNR低于20dB通常图像不可接受。
另外，由于数据集为彩色图，在图像压缩时会出现某些像素点明显失真的情况，表现为画面某些像素点出现明显错误的颜色，如下图：
这一现象可以由图像压缩中的近似误差来解释，由于RGB三种颜色的通道每一个像素点都是由一个0~255的数字表示的，在运算过程中每一步都会对256取模，这会导致某些像素点在0附近时会被近似到255，使得该通道原本无色的像素点变为纯色。
参考文献
[1]B. Qiu, V. Prinet, E. Perrier and O. Monga, “Multi-block PCA method for image change detection,” 12th International Conference on Image Analysis and Processing, 2003.Proceedings., 2003, pp. 385-390, doi: 10.1109/ICIAP.2003.1234080.
[2]S. T. Lim, D. F. W. Yap and N. A. Manap, “Medical image compression using block-based PCA algorithm,” 2014 International Conference on Computer, Communications, and Control Technology (I4CT), 2014, pp. 171-175, doi: 10.1109/I4CT.2014.6914169.
[3]D. Marvadi, C. Paunwala, M. Joshi and A. Vora, “Comparative analysis of 3D face recognition using 2D-PCA and 2D-LDA approaches,” 2015 5th Nirma University International Conference on Engineering (NUiCONE), 2015, pp. 1-5, doi: 10.1109/NUICONE.2015.7449603.
附代码main.py
# 对数据集应用三种PCA压缩方式并记录压缩效果、运行时间
import cv2
import numpy as np
from PIL import Image
from sklearn.metrics import mean_squared_error
from methods.pca_sklearn import *
from methods.pca_two_dimentions import *
from methods.pca_block_by_block import *
import time
def read_img(file_list, sample_num, size_n, size_m):
data = []
for filename in file_list:
img = cv2.imread(filename, 1)
# 预处理，调整尺寸为256px * 256px 插值填充
img = cv2.resize(img, (size_n, size_m), interpolation=cv2.INTER_AREA)
img_arr = np.asarray(img, dtype = int)
data.append(img_arr)
data = np.asarray(data, dtype = np.uint8)
return data
def save_images(pic_name, sample_num, compressed_data, size_n, size_m):
# 保存图像
for i in range(sample_num):
# 注意rgb顺序 实际上是bgr
img_b = compressed_data[i,:,:,0]
img_g = compressed_data[i,:,:,1]
img_r = compressed_data[i,:,:,2]
img = np.ndarray((size_n, size_m, 3), dtype = np.uint8)
img[:,:,0] = img_r
img[:,:,1] = img_g
img[:,:,2] = img_b
# print(img.shape)
image = Image.fromarray(img, 'RGB')
save_path = 'result/' + pic_name + '/' + pic_name + "%02d.tif" % i
image.save(save_path)
def calc_MSE_and_PSNR(data, compressed_data, sample_num, size_n, size_m):
# 对三个通道和100张图分别计算PSNR，取平均
MSE_SUM = np.mean((data - compressed_data) ** 2)
MSE = MSE_SUM
PSNR = 10 * np.log10(255 ** 2 / MSE)
return MSE, PSNR
if (__name__ == '__main__'):
pic_name = 'airplane'
sample_num = 100
# variance_ratio = 0.9
CR_list = [0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99]
size_n = 256
size_m = 256
file_list = []
for i in range(sample_num):
filename = 'data/' + pic_name + "/" + pic_name + "%02d.tif" % i
file_list.append(filename)
# 读取三个通道的数据(data[100,256,256,3])
data = read_img(file_list, sample_num, size_n, size_m)
f = open('result-%s.txt' % pic_name, 'w', encoding = 'utf-8')
f.write('Test started.\n')
f.write('datasets: %s\n' % pic_name)
f.write('method: traditional_pca\n')
print('method: traditional_pca')
for CR in CR_list:
start = time.time()
compressed_data, compression_ratio = compression_sklearn_pca_by_CR(data, CR)
# compressed_data, compression_ratio = two_dimention_pca_by_CR(data, CR)
# compressed_data, compression_ratio = compression_block_pca_by_CR(data, CR, 8)
MSE, PSNR = calc_MSE_and_PSNR(data, compressed_data, sample_num, size_n, size_m)
t = time.time() - start
f.write('CR: %.2f\t compression_ratio: %.2f\t MSE: %.2f\t PSNR: %.2f\t time: %.2fs\n' % (CR, compression_ratio, MSE, PSNR, t))
# save_images(pic_name, sample_num, compressed_data, size_n, size_m)
f.write('method: block_by_block_pca; block_size = 8\n')
print('method: block_by_block_pca; block_size = 8')
for CR in CR_list:
start = time.time()
# compressed_data, compression_ratio = compression_sklearn_pca_by_CR(data, CR)
# compressed_data, compression_ratio = two_dimention_pca_by_CR(data, CR)
compressed_data, compression_ratio = compression_block_pca_by_CR(data, CR, 8)
MSE, PSNR = calc_MSE_and_PSNR(data, compressed_data, sample_num, size_n, size_m)
t = time.time() - start
f.write('CR: %.2f\t compression_ratio: %.2f\t MSE: %.2f\t PSNR: %.2f\t time: %.2fs\n' % (CR, compression_ratio, MSE, PSNR, t))
f.write('method: two_dimentions_pca\n')
print('method: two_dimentions_pca')
for CR in CR_list:
start = time.time()
# compressed_data, compression_ratio = compression_sklearn_pca_by_CR(data, CR)
compressed_data, compression_ratio = two_dimention_pca_by_CR(data, CR)
# compressed_data, compression_ratio = compression_block_pca_by_CR(data, CR, 8)
MSE, PSNR = calc_MSE_and_PSNR(data, compressed_data, sample_num, size_n, size_m)
t = time.time() - start
f.write('CR: %.2f\t compression_ratio: %.2f\t MSE: %.2f\t PSNR: %.2f\t time: %.2fs\n' % (CR, compression_ratio, MSE, PSNR, t))
methods/pca_block_by_block.py
import numpy as np
def block_pca(data, block_size, variance_ratio):
data_mean = np.mean(data, axis = 1, keepdims=True)
data -= data_mean
C = data.T @ data / (100 - 1)
eig, V = np.linalg.eig(C)
V = V.T
idx = np.argsort(eig)[::-1]
# print(idx)
eig = eig[idx]
eig_sum = sum(eig)
V = V[idx]
# variance_ratio = 0.95
variance_sum = 0
n_components = 0
W = []
while (variance_sum / eig_sum < variance_ratio):
variance_sum += eig[n_components]
W.append(V[n_components])
n_components += 1
# print(n_components)
W = np.asarray(W)
W = W.T
z = W.T @ data.T
new_data = (W @ z).T
new_data += data_mean
ori_size = block_size * block_size * 100
com_size = n_components * 100 + n_components * block_size * block_size
return new_data.reshape(100, block_size, block_size), com_size / ori_size
def block_pca_by_n_c(data, block_size, n_c):
data_mean = np.mean(data, axis = 1, keepdims=True)
data -= data_mean
C = data.T @ data / (100 - 1)
eig, V = np.linalg.eig(C)
V = V.T
idx = np.argsort(eig)[::-1]
# print(idx)
eig = eig[idx]
eig_sum = sum(eig)
V = V[idx]
# variance_ratio = 0.95
variance_sum = 0
n_components = 0
W = []
while (n_components < n_c):
variance_sum += eig[n_components]
W.append(V[n_components])
n_components += 1
# print(n_components)
W = np.asarray(W)
W = W.T
z = W.T @ data.T
new_data = (W @ z).T
new_data += data_mean
ori_size = block_size * block_size * 100
com_size = n_components * 100 + n_components * block_size * block_size
# CR = (n_components * 100 + n_components * block_size * block_size) / (block_size * block_size * 100)
# -> n_c = (CR * bs * bs * 100) / (100 + bs * bs)
return new_data.reshape(100, block_size, block_size), com_size / ori_size
def compression_block_pca(data, variance_ratio, block_size): # data[100, 256, 256, 3]
print('compression_block_pca:')
# block_size = 8
n_samples = data.shape[0]
size_n = data.shape[1]
size_m = data.shape[2]
n_channel = data.shape[3]
n_row = size_n // block_size
n_col = size_m // block_size
data_total = np.ndarray((n_row, n_col, 100, block_size * block_size, 3))
# print(data.shape)
for i in range(n_channel):
for r in range(n_row):
for c in range(n_col):
for j in range(n_samples):
ur = r * block_size
dr = (r + 1) * block_size
lc = c * block_size
rc = (c + 1) * block_size
data_total[r, c, j, :, i] = np.asarray(data[j, ur:dr, lc:rc, i]).reshape(1, block_size * block_size)
# print(data_total.shape)
compressed_data = np.ndarray((100, 256, 256, 3))
ratio_sum = 0
cnt = 0
for i in range(n_channel):
for r in range(n_row):
for c in range(n_col):
data_block = data_total[r, c, :, :, i]
ur = r * block_size
dr = (r + 1) * block_size
lc = c * block_size
rc = (c + 1) * block_size
# print(data_block.shape)
compressed_data[:, ur:dr, lc:rc, i], ratio = block_pca(data_block, block_size, variance_ratio)
ratio_sum += ratio
# print(cnt)
cnt += 1
return compressed_data, ratio_sum / (n_channel * n_row * n_col)
def compression_block_pca_by_CR(data, CR, block_size): # data[100, 256, 256, 3]
print('compression_block_pca:')
# block_size = 16
n_samples = data.shape[0]
size_n = data.shape[1]
size_m = data.shape[2]
n_channel = data.shape[3]
n_row = size_n // block_size
n_col = size_m // block_size
data_total = np.ndarray((n_row, n_col, 100, block_size * block_size, 3))
n_c = int((CR * block_size * block_size * 100) / (100 + block_size * block_size)) + 1
print(CR, n_c)
# print(data.shape)
for i in range(n_channel):
for r in range(n_row):
for c in range(n_col):
for j in range(n_samples):
ur = r * block_size
dr = (r + 1) * block_size
lc = c * block_size
rc = (c + 1) * block_size
data_total[r, c, j, :, i] = np.asarray(data[j, ur:dr, lc:rc, i]).reshape(1, block_size * block_size)
# print(data_total.shape)
compressed_data = np.ndarray((100, 256, 256, 3))
ratio_sum = 0
cnt = 0
for i in range(n_channel):
for r in range(n_row):
for c in range(n_col):
data_block = data_total[r, c, :, :, i]
ur = r * block_size
dr = (r + 1) * block_size
lc = c * block_size
rc = (c + 1) * block_size
# print(data_block.shape)
compressed_data[:, ur:dr, lc:rc, i], ratio = block_pca_by_n_c(data_block, block_size, n_c)
ratio_sum += ratio
# print(cnt)
cnt += 1
return compressed_data, ratio_sum / (n_channel * n_row * n_col)
methods/pca_two_dimentions.py
import numpy as np
def two_dimention_pca(data, variance_ratio): # data[100, 256, 256, 3]
print('two_dimention_pca:')
n_samples = data.shape[0]
size_n = data.shape[1]
size_m = data.shape[2]
n_channel = data.shape[3]
compressed_data = np.ndarray((n_samples, size_n, size_m, n_channel))
ratio = 0
for i in range(n_channel):
average = np.zeros((size_n, size_m))
for j in range(n_samples):
average += data[j, :, :, i]
average /= n_samples
G_t = np.zeros((size_m, size_m))
for j in range(n_samples):
img = data[j, :, :, i]
temp = img - average
G_t = G_t + np.dot(temp.T, temp) / n_samples
eig, V = np.linalg.eig(G_t)
idx = np.argsort(eig)
idx = idx[::-1]
eig_sum = sum(eig)
eig = eig[idx]
V = V.T
V = V[idx]
variance_sum = 0
n_components = 0
W = []
while (variance_sum / eig_sum < variance_ratio):
variance_sum += eig[n_components]
W.append(V[n_components])
n_components += 1
# print(n_components)
W = np.asarray(W)
W = W.T
# print(W.shape)
newx = (data[:,:,:,i] - average) @ W
X = newx @ W.T + average
# print(X.shape)
ori_size = n_samples * size_n * size_m
com_size = n_components * n_samples * size_m + n_components * size_m
ratio = com_size / ori_size
compressed_data[:,:,:,i] = X
ratio /= n_channel
return compressed_data, ratio
def two_dimention_pca_by_CR(data, CR): # data[100, 256, 256, 3]
print('two_dimention_pca:')
n_samples = data.shape[0]
size_n = data.shape[1]
size_m = data.shape[2]
n_channel = data.shape[3]
compressed_data = np.ndarray((n_samples, size_n, size_m, n_channel))
n_c = int(CR * n_samples * size_n * size_m / (n_samples * size_m + size_m)) + 1
print(CR, n_c)
ratio = 0
for i in range(n_channel):
average = np.zeros((size_n, size_m))
for j in range(n_samples):
average += data[j, :, :, i]
average /= n_samples
G_t = np.zeros((size_m, size_m))
for j in range(n_samples):
img = data[j, :, :, i]
temp = img - average
G_t = G_t + np.dot(temp.T, temp) / n_samples
eig, V = np.linalg.eig(G_t)
idx = np.argsort(eig)
idx = idx[::-1]
eig_sum = sum(eig)
eig = eig[idx]
V = V.T
V = V[idx]
variance_sum = 0
n_components = 0
W = []
while (n_components < n_c):
variance_sum += eig[n_components]
W.append(V[n_components])
n_components += 1
# print(n_components)
W = np.asarray(W)
W = W.T
# print(W.shape)
newx = (data[:,:,:,i] - average) @ W
X = newx @ W.T + average
# print(X.shape)
ori_size = n_samples * size_n * size_m
com_size = n_components * n_samples * size_m + n_components * size_m
# CR = (n_components * n_samples * size_m + n_components * size_m) / (n_samples * size_n * size_m)
# -> n_components = CR * n_samples * size_n * size_m / (n_samples * size_m + size_m)
ratio += com_size / ori_size
compressed_data[:,:,:,i] = X
ratio /= n_channel
return compressed_data, ratio