转载

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

上一个笔记主要是讲了PCA的原理，并给出了二维图像降一维的示例代码。但还遗留了以下几个问题：

在计算协方差和特征向量的方法上，书上使用的是一种被作者称为 compact trick 的技巧，以及 奇异值分解（SVD） ，这些都是什么东西呢？
如何把PCA运用在多张图片上？

所以，我们需要进一步的了解，同时，为示例对多张图片进行PCA，我选了一个跟书相似但更有趣的例子来做——人脸识别。

人脸识别与特征脸

一个 特征脸（Eigenface，也叫标准脸） 其实就是从一组人脸图像应用PCA获得的主成分特征向量之一，下面我们能验证，每个这样的特征向量变换为二维图像后看起来有点像人脸，所以才被称为 特征脸 ，计算特征脸是进行人脸识别的首要步骤，其计算过程其实就是PCA。

特征脸计算步骤

准备一组（假设10张）具有相同分辨率（假设:100 × 100）的人脸图像，把每张图像打平（numpy.flatten）成一个向量，即所有像素按行串联起来, 每个图像被当作是10000维空间的一个点。再把所有打平的图像存储在 10000 × 10 矩阵X中，矩阵的每一列就是一张图片，每个维为一行，共10000维
对X的每一维（行）求均值，得到一个10000 × 1的向量，称为 平均脸 （因为把它变换为二维图像看起来像人脸），然后将X减去平均脸，即零均值化
计算X的协方差矩阵 C=XX^(X^表示X的转置)
计算C的特征值和特征向量，这组特征向量就是一组特征脸。在实际应用中，我们只需要保留最主要的一部分特征向量做为特征脸即可。

有了上个笔记的基础知识后，上面计算过程不难理解。但在实现代码之前，我们先来看看上面提到的计算X的协方差矩阵 C=XX^ 引发的一个问题。

协方差计算技巧

对于上面举例的矩阵X，它有10000行（维），它的协方差矩阵将达到 10000 × 10000，有10000个特征向量，这个计算量是很大的，消耗内存大，我尝试过不可行。

这个数学问题终究还是被数学解决了，解决的方法可以见维基百科 特征脸 内容描述。简言之，就是把本来计算 XX^ 的协方差矩阵（设为C）变换为计算 X^X 的协方差矩阵（C'），后者的结果是 10 × 10（10为样本数量），很快就可以算出来。当然，通过这个变换分别算出来的特征向量不是等价的，也需要变换一下：设E为从C算出来的特征向量矩阵，E'为从C'算出来的特征向量矩阵，则 E = XE' ，最后再把E 归一化 。

这个技巧就是书上PCA示例使用的，被称为 compact trick 的方法。

但要看明白书上的示例代码，还要搞清一点：

对原始图像数据集矩阵的组织方式，我们用行表示维度，列表示样本，而书上和《Guide to face recognition with Python》（见底部参考链接）使用的是行表示样本，列表示维度。就是因为这两种组织方式的不同，导致了PCA算法的代码看起来有些不同。这一点很容易让人困惑，所以写到这里，我应该特别的强调一下。

我之所以在上个笔记，包括上面对特征脸的计算步骤描述，都认定以行表示维度，列表示样本的方式，是为了与数学原理的详解保持一致（注：下面的代码示例还是使用这种方式），当我们明白了整个原理之后，我们便知道使用这两种矩阵表达方式都可以，两者实现的PCA代码差别也很小，下面会讲到。

准备人脸图像样本

网上有不少用于研究的人脸数据库可以下载，我在参考链接给出了常被使用的一个。下载解压后在目录orl_faces下包含命名为s1,..,s40共40个文件夹，每个文件夹对应一个人，其中存储10张脸照，所有脸照都是92 × 112的灰度图，我把部分照片贴出来：

接下来，我们按照 特征脸计算步骤 中的第1点所述，把这400张图像组成矩阵（图像组织不分先后），代码：

def getimpaths(datapath):  paths = []  for dir in os.listdir(datapath):   try:    for filename in os.listdir(os.path.join(datapath, dir)):     paths.append(os.path.join(datapath, dir, filename))   except:    pass  return paths impaths = getimpaths('./orl_faces') m,n = np.array(Image.open(impaths[0])).shape[0:2] #图片的分辨率，下面会用到 X = np.mat([ np.array(Image.open(impath)).flatten() for impath in impaths ]).T print 'X.shape=',X.shape  #X.shape= (10304, 400)

我们把每个图像都打平成行向量，把所有图像从上到下逐行排列，最后转置一下，便得到一个10304 × 400 的矩阵，其中10304 = 92 × 112

实现PCA函数求解特征脸

PCA函数

我尽量使用与书上相同的变量命名，方便对比：

def pca(X):     dim, num_data = X.shape  #dim: 维数， num_data: 样本数          mean_X = X.mean(axis=1) #求出平均脸，axis=1表示计算每行（维）均值，结果为列向量     X = X - mean_X          #零均值化          M = np.dot(X.T, X)         #使用compact trick技巧，计算协方差     e,EV = np.linalg.eigh(M) #求出特征值和特征向量     print 'e=',e.shape,e     print 'EV=',EV.shape,EV          tmp = np.dot(X, EV).T      #因上面使用了compact trick，所以需要变换     print 'tmp=',tmp.shape,tmp     V = tmp[::-1]              #将tmp倒序，特征值大的对应的特征向量排前面，方便我们挑选前N个作为主成分     print 'V=',V.shape,V          for i in range(num_data):         V[:,i] /= np.linalg.norm ( EV[:,i]) #因上面使用了compact trick，所以需要将V归一化          return V,EV,mean_X

执行PCA并画图对上面得到的X矩阵调用pca函数，并画出平均脸和部分特征脸：

V,EV,immean = pca(X)  #画图 plt.gray() plt.subplot(2,4,1)  #2行4列表格，第一格显示`平均脸` plt.imshow(immean.reshape(m,n))  #以下选前面7个特征脸，按顺序分别显示到其余7个格子 for i in range(7):     plt.subplot(2,4,i+2)     plt.imshow(V[i].reshape(m,n))      plt.show()

显示效果图如下：

希望不会被这些特征脸吓到：）

这些所谓的脸事实上是特征向量，只不过维数与原始图像一致，因此可以被变换成图像显示出来，不同的特征脸代表了与均值图像差别的不同方向。

当然，我们求特征脸，并不是为了显示他们，而是保留部分特征脸来获得大多数脸的近似组合。因此，人脸便可通过一系列向量而不是原始数字图像进行保存，节省了很多存储空间，也便于后续的识别计算。

与书上pca的实现对比

上面我给出的pca函数代码，是按照我们一路学习PCA的思路写出来的，虽然跟书上pca实现很接近，但还是有几个点值得分析：

如上提到，我们对X矩阵的组织是以行为维、列为样本的方式，即一个列对应一张打平的图像，而书上的例子是以行为样本、列为维的方式，每一行对应一张打平的图像，而且参考链接里的例子也都是以后者进行组织的，但没关系，我们只需要对上面的代码作一点修改即可：

def pca_book(X):     num_data, dim = X.shape     #注意：这里行为样本数，列为维          mean_X = X.mean(axis=0) #注意：axis=0表示计算每列（维）均值，结果为行向量     X = X - mean_X                      #M = np.dot(X.T, X)            #把X转置后代入，得到     M = np.dot(X, X.T)            #跟书上一样          e,EV = np.linalg.eigh(M)    #求出特征值和特征向量          #tmp = np.dot(X, EV).T        #把X转置后代入，得到     tmp = np.dot(X.T, EV).T     #跟书上一样      V = tmp[::-1]                          #以下是对V归一化处理，先省略，下面讲

所以我们看到，其实算法的本质是一样的，只不过要注意的地方是维数和样本数反过来了，另外，对X的运算换成X的转置即可。类推的，如果X使用我们的上面的组织方式，调用pca_book函数的代码为 V,EV,immean = pca_book(X.T)

归一化 算法不同。因为使用书上的方法，在对特征值求平方根( np.sqrt(e) )的时候会产生一个错误：有个别数无法计算平方根，所以我从参考链接例子获得的方法是使用 np.linalg.norm 这个函数
书上的pca算法多了一个判断分支，当 dim <= num_data 即维数少于样本数的时候直接使用 SVD（奇异值分解） 算法，显然在一般的人脸识别的例子里，不会被用到，因为单张92 × 112图像打平后维数为10304，而样本数为400，远远低于维数。

归一化

待写。

SVD（奇异值分解）

回顾一下上篇笔记举的PCA应用例子：把一张二维图像变换成一维，维数为2，而样本数为所有黑点像素，很大，对于这个例子，直接使用SVD是比较合适的，这样PCA函数将变得很简洁。