Programming Computer Vision with Python （学习笔记三）

概要

原书对于PCA的讲解只有一小节，一笔带过的感受，但我发现PCA是一个很重要的基础知识点，在机器机视觉、人脸识别以及一些高级图像处理技术时都被常常用到，因此本人自行对PCA进行了更深刻的学习。

PCA是什么

PCA（Principal Component Analysis，主成分分析或主元分析）是一种算法，PCA的结果是用尽量少的特征数据来表达最多的原始图像的本质结构特征。即便它会丢失一部分原始图像的特征表达，但它仍然是颇有用的处理技巧，也很经常使用，特别在计算机视觉和人脸识别方面。

假设咱们有一个二维数据集，它在平面上的分布以下图：

若是咱们想要用一个一维向量来表达此数据集，就会丢失一部分此数据集的信息，但咱们的目标是让求得的这个一维向量能够尽量多地保留这个数据集的特征信息，那么这个求解过程就是PCA。

经过PCA咱们能够找到若干个1维向量，如图：

直观上看出，向量u1是数据集变化的主方向，而u2是次方向，u1比u2保留了更多的数据集的结构特征，因此咱们选择u1做为主成分，并把原数据集投影到u1上就能够得出对原数据集的一维近似重构：

以上只是一种直观的示例，把二维数据降为用1维来表示，固然，PCA一般是应用在高维数据集上。

PCA解决什么问题

假设咱们有10张100 × 100像素的灰度人脸图，咱们目标是要计算这10张图的主成分来做为人脸特征，这样就能够基于这个‘特征脸’进行人脸匹配和识别。但即便一个100 × 100像素的灰度图像维度就达到10,000维，10张图像的线性表示能够达到100,000维，如此高维的数据带来几个问题：

• 对高维数据集进行分析处理的计算量是巨大的，消耗资源太大，时间太长

• 高维数据包含了大量冗余和噪声数据，会下降图像识别率

因此一般对于高维数据集，首先须要对其进行降维运算，以低维向量表达原数据集最多最主要的结构特征。从而将高维图像识别问题转化为低维特征向量的识别问题，大大下降了计算复杂度，同时也减小了冗余信息所形成的识别偏差。PCA其实就是最经常使用的一种降维算法。

PAC也可用于高维数据压缩、高维数据可视化（转二维或三维后就能够画图显示）等方面，也是其它不少图像处理算法的预处理步骤。

PCA的计算

关于PCA，网上一搜仍是很多的，但我仔细看了几篇文章以后，发现这些文章讲的跟书上讲的有些地方不一致，甚至连计算时的公式都不同，这让我产生了不少困惑。因此我想最重要的仍是要理解PCA的数学原理，数学原理才是根，掌握了数学原理，再来写代码。刚好找到一篇文章专门讲PCA数学原理，做者的数学功底和逻辑表达能力很是棒，让我很容易看明白。另外，我也找到一篇老外写的文章（见底部参考文章），这两篇文章对PCA的计算描述是一致的，因此我决定在这两篇文章的基础上，结合书上的示例代码进行学习和验证。

本文不是要要把PCA的数学原理及推导写出来，而是经过理解PCA的数学原理，总结PCA的计算步骤，由于计算步骤是代码实现的必备知识。

PCA的计算过程涉及到几个很重要的数学知识点：

• 零均值化

• 矩阵的转置及乘法

• 协方差与协方差矩阵

• 特征值及特征向量

如今来看PCA的计算步骤：
1）将原始数据按列组成d行n列矩阵X
重要说明：d对应的就是数据的字段（或叫变量、特征、维，下称’维‘），而n表示n条记录（或叫样本、观察值，下称’样本‘），即每1列对应1个样本，之因此行和列这样安排，是为了与数学公式保持一致，不少文章对这一点都没有明确的说明，致使计算的方式各有不一样，让人产生没必要要的困惑
2）将X的每一个维（行）进行零均值化，即将行的每一个元素减去这一行的均值
3）求出X的协方差矩阵C，即 X 乘 X的转置
4）求出C全部的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵E
6）Y=EX即为降维到k维后的数据

下面用一个例子来验证一下这个计算过程。

黑白图像的PCA实现

书中的例子是用PCA计算特征脸（人脸识别中的一步），它应用在多张图片上面。为直观起见，我用另一个例子——对单张黑白图像进行PCA，至关于把二维图像降为一维。

把图像转成二维矩阵
这张图像是原书封面：

下面代码是将图中‘怪鱼’部分截取出来，并转成黑白图像显示：

from PIL import Image
pim = Image.open('cover.png').crop((110,360,460,675)).convert('1')
pim.show()

效果如图：

之因此截取这部分的图片，是由于咱们大概能猜到这幅图像降到一维后，其一维表示的向量应该跟怪鱼的方向大概一致。

使用黑白图像是由于黑点才是咱们关心的数据，由于是这些黑点描绘了图像，每一个黑点有惟一肯定的行和列位置，对应平面上的（x,y）坐标，因而咱们就能够获得此图像的 2乘n 矩阵表示：第一行表示x维，第二行表示y维，每一列表示一个点。参考代码：

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

im = np.array(Image.open('cover.png').crop((110,360,460,675)).resize((256,230)).convert('L'))
n,m = im.shape[0:2]
points = []
for i in range(n):
    for j in range(m):
        if im[i,j] < 128.0:  #把小于128的灰度值看成黑点取出来
            points.append([float(j), float(n) - float(i)]) #坐标转换一下
    
im_X = np.mat(points).T; #转置以后，行表示维度（x和y），每列表示一个点（样本）
print 'im_X=',im_X,'shape=',im_X.shape

如今，咱们按上面说明的计算步骤来实现PCA：

def pca(X, k=1): #降为k维
  d,n = X.shape
  mean_X = np.mean(X, axis=1) #axis为0表示计算每列的均值，为1表示计算每行均值
  print 'mean_X=',mean_X
  X = X - mean_X
  #计算不一样维度间的协方差，而不是样本间的协方差，方法1：
  #C = np.cov(X, rowvar=1) #计算协方差，rowvar为0则X的行表示样本，列表示特征/维度
  #方法2：
  C = np.dot(X, X.T)
  e,EV = np.linalg.eig(np.mat(C)) #求协方差的特征值和特征向量
  print 'C=',C
  print 'e=',e
  print 'EV=',EV
  e_idx = np.argsort(-e)[:k] #获取前k个最大的特征值对应的下标（注：这里使用对负e排序的技巧，反而让本来最大的排在前面）
  EV_main = EV[:,e_idx]   #获取特征值（下标）对应的特征向量，做为主成分
  print 'e_idx=',e_idx,'EV_main=',EV_main      
  low_X = np.dot(EV_main.T, X)    #这就是咱们要的原始数据集在主成分上的投影结果
  return low_X, EV_main, mean_X

OK,如今咱们调用此PCA函数，并把原图像和投影到一维向量后的结果也描绘出来：

low_X, EV_main, mean_X = pca(im_X)
print "low_X=",low_X
print "EV_main=",EV_main
recon_X = np.dot(EV_main, low_X) + mean_X  #把投影结果重构为二维表示，以即可以画出来直观的看到
print "recon_X.shape=",recon_X.shape

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(im_X[0].A[0], im_X[1].A[0],s=1,alpha=0.5)
ax.scatter(recon_X[0].A[0], recon_X[1].A[0],marker='o',s=100,c='blue',edgecolors='white')
plt.show()

画散点图函数pyplot.scatter说明：

matplotlib.pyplot.scatter(x, y, ...)
x: 数组，样本点在X轴上的坐标集
y: 数组，样本点在Y轴上的坐标集
s: 表示画出来的点的缩放大小
c: 表示画出来的点（小圆圈）的内部颜色
edgecolors: 表示小圆圈的边缘颜色

运行以上代码打印：

im_X= [[  23.   24.   25. ...,  215.  216.  217.]
 [ 230.  230.  230. ...,    5.    5.    5.]] shape= (2, 19124)
mean_X= [[ 133.8574566 ]
 [ 123.75941226]]
C= [[ 2951.65745054 -1202.25277635]
 [-1202.25277635  3142.71830026]]
e= [ 1841.14567037  4253.23008043]
EV= [[-0.73457806  0.67852419]
 [-0.67852419 -0.73457806]]
e_idx= [1] EV_main= [[ 0.67852419]
 [-0.73457806]]
low_X= [[-153.26147057 -152.58294638 -151.90442219 ...,  142.29523704
   142.97376123  143.65228541]]
EV_main= [[ 0.67852419]
 [-0.73457806]]
recon_X.shape= (2, 19124)

并显示以下图像：

从图中看出，向量的方向跟位置与咱们目测的比较一致。向量上的大量蓝色圆点（白色边缘）表示二维数据在其上的投影。

小结

以上实现的PCA算法，跟我参考的两篇文章所讲的原理一致。但跟书中的PCA计算方法有必定的不一样，但也由于使用的例子不同，对原始数据集的定义不同致使的，因为避免文章过长，这些放在后面再讲。
至此，经过对PCA的计算过程的学习，了解了一些线性代数知识、numpy和pyplot模块的一些接口的用法，接下来我打算作点更有兴趣的事情——就是使用PCA来实现人脸识别。

参考连接：
PCA的数学原理
A Tutorial on Principal Component Analysis
NumPy函数索引