PCA 实践 利用 PCA 算法对人脸数据集内所有人进行降维和特征提取 PCA原理解析+代码
PCA 实践 利用 PCA 算法对人脸数据集内所有人进行降维和特征提取 PCA原理解析+代码
数据集
实验所用到的数据集在下面的链接中， 这些数据是来自剑桥大学提供的 AT&T 人脸数据 集,有 40 个人的人脸图像, 每个人有 10 张不同光照和姿态的照片。
样例：
地址： http://www.cl.cam.ac.uk/Research/DTG/attarchive/pub/data/att_faces.tar.Z
实验内容
加载数据集，利用 PCA 算法对数据集内所有人进行降维和特征提取，然后将得到的主 成分特征向量还原成图像进行观察。这里可以尝试采用不同的降维维度 K 进行操作，分别观 察不同 K 下的特征图像。
实验拓展
尝试对刚降维的特征图像进行 PCA 逆变换，观察变换前后的图像差异
- 实验步骤与内容：
- 分析数据集
数据集中包含了40个文件夹，对应了40个志愿者人脸，每个文件夹中有一个志愿者的人脸的多种状态。
- 分析PCA原理：
- 假设在Rn空间中有m个点, 我们希望对这些点进行有损压缩, 使数据的维度从Rn变为Rl, 其中严格的有l<n. 这时候主成分分析法(PCA)便可以实现我们的要求.对于每个点xi∈Rn, 使得其被投影为ci∈Rl, 用函数来表示下列编码过程即是: f(x)=DTx=c。 同时, 我们也希望找到一个解码函数, 使得 g(f(x))=g(c)=Dc=DDTx≈x′其中D是一个列向量彼此正交的Rn∗lRn∗l矩阵. PCA有两种推导过程, 但它们的结论是一样的。
- PCA的两种推导过程：1.最大方差理论；2.最小误差理论；
- PCA算法步骤：
- 关于方差和协方差
- 协方差矩阵
- 读取并处理数据集
此处使用相对目录，依次访问每个子文件夹，通过在中间过程中输出图像的维度，我们观察到图像是112x92x3即，长宽分别为112,93，通道数为3。首先，我们先对图像进行灰度化处理，然后将图像维度展开成为一列，存入列表。
- 图像显示测试
图像的height，width=112,92，X代表我们的数据集，是一个400x10304的向量；之所以是400x10304，是因为总共400幅图像，每个图像是一个112x92的image，将图像展开成一个一维向量，维度为1x10304，将400幅图像组合得到一个400x10304的矩阵。其中Y是数据的标签。
此处，我们展示试着显示一组图像，在文件夹中选取每个子文件夹下的第一个图像拼合显示如下：
因为窗口大小限制，此处并没有显示完全40副图像。
- 计算平均脸
其中，对10304个列向量做平均，即对400个人脸图像的每一个对应像素值做平均，计算均值，得到一个1x10304的向量。并将该向量重新转化为和原图像同样大小的维度，然后显示出来，如下：
- 数据中心化，并进行奇异值分解
注意，计算协方差矩阵之前，需要先将每个变量减去均值，得到中心化的向量，然后通过使用sklearn中的奇异值分解函数，得到前K个特征向量。
- 之所以使用奇异值分解而不是协方差矩阵的特征值分解，是因为如下原因：
不需要计算协方差矩阵S，同时，在数值上更加精确，因为在计算机存储中，可能会产生累计误差。
- K值的确定：
我们设定阈值为0.95，计算保留的主成分数量。计算代码如下：
计算得到的K值为111，也就是说，前111个奇异值得到的结果已经能代表90%的情况了，为方便处理，我们此处选取K=100.通过奇异值分解得到Kx10304的变换矩阵，我们可以得到如下十张主成分图像：
注意：如果使用OpenCV中的imshow函数进行显示，需要先将数据转换为np.unit8格式，同时范围规范化到0-255的范围（此处使用最大最小值归一化），而使用matplotlib则不用进行如此操作。
规范化及使用OpenCV显示相关代码
使用matplotlib显示相关代码：
- 图像还原：
当完成主成分提取之后，我们对图像进行还原：
其中，U是变换矩阵，还原方式如下：
并输出30张还原之后的图像；如下：
我们可以看到，以上30张图像分别对应了原始数据集中的3个人的图像，相比于原来的30张，如今还原得到的图像丢失了一部分信息，不过不是非常明显。
源代码
1. import cv2
2. import tensorflow as tf
3. import numpy as np
4. import matplotlib.pyplot as plt
5. from sklearn.decomposition import TruncatedSVD
6.
7. FACE_PATH = "orl_faces"
8. PERSON_NUM = 40
9. PERSON_FACE_NUM = 10
10. K = 10 # Number of principle components
11.
12. raw_img = []
13. data_set = []
14. data_set_label = []
15.
16.
17. def read_data():
18. height = 0
19. width = 0
20. for i in range(1, PERSON_NUM + 1):
21. person_path = FACE_PATH + '/s' + str(i)
22.
23. for j in range(1, PERSON_FACE_NUM + 1):
24. img = cv2.imread(person_path + '/' + str(j) + '.pgm')
25. #print(img.shape) 112,92,3
26. if j == 1:
27. raw_img.append(img)
28. img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
29. height, width = img_gray.shape
30. #print(height,width)
31. img_col = img_gray.reshape(height * width)
32. data_set.append(img_col)
33. data_set_label.append(i)
34.
35. return height, width
36.
37.
38. # Import Data
39. height, width = read_data()
40. print(height,width)
41. X = np.array(data_set)
42. #print(X.shape) # 400*10304
43. Y = np.array(data_set_label)
44. #print(Y.shape)
45. n_sample, n_feature = X.shape
46.
47. # Print some samples
48. raw_img = np.hstack(raw_img)
49. cv2.namedWindow("Image")
50. cv2.imshow('Image', raw_img)
51. cv2.waitKey(1)
52. cv2.destroyAllWindows()
53.
54. # 计算均值 A preview of average_face
55. average_face = np.mean(X, axis=0) #对10304个列向量做平均
56. # for i in average_face:
57. # print(i)
58. #num+=1
59. fig = plt.figure()
60. plt.imshow(average_face.reshape((height, width)), cmap=plt.cm.gray)
61. plt.title("Average Face", size=12)
62. plt.xticks(())
63. plt.yticks(())#不显示横纵坐标
64. #plt.savefig("average_face.png")
65. plt.show()
66.
67. # Calculate
68. equalization_X = X - average_face #差值400*10304
69. #计算协方差矩阵
70.
71. covariance_X = np.cov(equalization_X.transpose())
72. print(covariance_X.shape)
73. print(equalization_X.shape)
74. #特征值分解 使用sklearn中的SVD函数分解
75. svd = TruncatedSVD(n_components=K, random_state=44)
76. svd.fit(equalization_X)
77.
78. # total_variance, n_components = 0.0, 0
79. # for variance in svd.explained_variance_ratio_:
80. # total_variance += variance
81. # n_components += 1
82. # if total_variance > 0.95: break
83.
84. #print(n_components) 9
85. #print(svd.components_.shape) 10*10304
86. print(svd.components_.shape)
87. topk=np.zeros(svd.components_.shape)
88. # for i in range(K):
89. # #print(svd.components_[i])
90. # topk[i]=255.0*(svd.components_[i]-np.min(svd.components_[i]))/\
91. # (np.max(svd.components_[i])-np.min(svd.components_[i]))
92. #print(topk[i])
93. # img=topk[1].reshape(height, width)
94. # img=img.astype(np.uint8)
95. # #print(img.shape)
96. # cv2.imshow('1',img)
97. # cv2.waitKey(0)
98. # cv2.destroyAllWindows()
99. #np.matmul(equalization_X,)
100. # topk=svd.components_+average_face
101.
102.
103. # 打印输出前K个主成分脸
104.
105. plt.figure()
106. for i in range(1, K + 1):
107. plt.subplot(2, K/2, i)
108. #print(svd.components_[i-1].shape)
109. plt.imshow((svd.components_[i - 1]).reshape(height, width), cmap=plt.cm.gray)
110. #print(svd.components_[i - 1].reshape(height, width))
111. plt.xticks(())
112. plt.yticks(())
113. plt.show()
114.
115. topk=average_face+np.matmul(equalization_X,np.matmul(svd.components_.T,svd.components_))
116. print(topk.shape)
117. K=30
118. plt.figure()
119. for i in range(1, K + 1):
120. plt.subplot(5, K/5, i)
121. #print(svd.components_[i-1].shape)
122. plt.imshow((topk[i - 1]).reshape(height, width), cmap=plt.cm.gray)
123. #print(svd.components_[i - 1].reshape(height, width))
124. plt.xticks(())
125. plt.yticks(())
126. plt.show() 如不知道如何去掉行号，参考个人博文中的解决方案。