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

摘要：从而将高维图像识别问题转化为低维特征向量的识别问题，大大降低了计算复杂度，同时也减少了冗余信息所造成的识别误差。向量上的大量蓝色圆点白色边缘表示二维数据在其上的投影。

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

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

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

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

通过PCA我们可以找到若干个1维向量，如图：

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

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

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

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

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

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

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

关于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，相当于把二维图像降为一维。

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

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

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函数索引