实验所用到的数据集在下面的链接中, 这些数据是来自剑桥大学提供的 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个志愿者人脸,每个文件夹中有一个志愿者的人脸的多种状态。
此处使用相对目录,依次访问每个子文件夹,通过在中间过程中输出图像的维度,我们观察到图像是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,同时,在数值上更加精确,因为在计算机存储中,可能会产生累计误差。
我们设定阈值为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()
如不知道如何去掉行号,参考个人博文中的解决方案。