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函数
我尽量使用与书上相同的变量命名，方便对比：

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(EV.shape[1]):
        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))的时候会产生一个错误（负数不能开平方根），所以我这里使用的归一化方法是从《Guide to face recognition with Python》抄来的。

书上的pca算法多了一个判断分支，当dim <= num_data即维数少于样本数的时候直接使用SVD（奇异值分解）算法，显然在一般的人脸识别的例子里，不会被用到，因为单张92 × 112图像打平后维数为10304，而样本数为400，远远低于维数。

归一化
原先以为归一化是一种比较简单的运算，一了解才发现原来是一种不简单的思想，在机器学习中常被使用，看回上面的代码：

for i in range(EV.shape[1]):
        V[:,i] /= np.linalg.norm ( EV[:,i])

首先读者得自行了解范数(norm)的概念, 范数（norm）还分L0、L1、L2等好几种，而函数np.linalg.norm就是用来求矩阵或向量的各种范数，默认就是求L2范数，具体可查阅《linalg.norm API说明》
上面代码的作用就是对V每一列归一化到单位L2范数。
而书上使用的归一化方法是：

S = sqrt(e)[::-1] #计算e的平方根再对结果倒序排列
for i in range(V.shape[1]):
V[:,i] /= S

我在网上找到了关于compat trick后如何对求得的向量归一化的数学推荐，截图如下：

这跟左奇异值有关，属于SVD中的内容，有兴趣的话自行研究。
当我使用这种方法实现时，程序运行出现错误：FloatingPointError: invalid value encountered in sqrt，发现是对负数开平方根产生了错误，也就是说对协方差矩阵求得的特征值中包含了负数。那么，如果要使用这种归一化方法的话，是否只要排除掉负的特征值及其对应的特征向量就可以了？

SVD（奇异值分解）
我们的代码示例的PCA方法使用的是特征分解，线性代数中，特征分解（Eigendecomposition），又称谱分解（Spectral decomposition）,是将矩阵分解为由其特征值和特征向量表示的矩阵之积的方法，但需要注意只有对可对角化矩阵才可以施以特征分解。
而SVD（singular value decomposition）能夠用于任意m乘n矩阵的分解，故SVD适用范围更广。
但是，如果矩阵维数很大，如我们之前所举例10000维的时候，计算SVD也是很慢的，所以我们看到书上的例子增加了一个分支判断，当维度 < 样本数的时候，才使用SVD，否则使用compat trick方法的PCA。
回顾一下上篇笔记举的PCA应用例子：把一张二维图像变换成一维，维数为2，对于这个例子，直接使用SVD是比较合适的，这样PCA函数将变得很简洁。

这里不会详细地讨论如何实现一个好的人脸识别算法，而是为了示例PCA的运用，所以只是简单的介绍一下。
上面我们已经得出了400张人脸的特征脸（特征向量），首先第一个问题，我们得选择多少个特征向量作为主成分？
特征值贡献率
假设我们选择k个特征向量，其对应的特征值之和与所有特征值之和的比值，就是这k个特征值的贡献率。所以主成分的选择问题就转化为选择k个特征向量，使得特征值的贡献率大于等于某个值（如90%）即可。我们把选定的k个特征向量组成的矩阵设为W。

识别步骤
第一步：将每个人脸样本图像减去平均图像后，投影到主成分上

W = EV[:,k]       #假设k已经根据特征值贡献率算出来了
projections = []    #存放每个人脸样本的投影
for xi in X.T:    #X为我们之前组织的所有人脸样本的 10000  × 400矩阵
    projections.append(np.dot(W.T, xi - mean_X)) #mean_X为之前我们求得的平均图像

第二步：设要识别的图像为D，将D也投影到主成分上得到Q，然后计算Q与各个样本人脸投影的欧几里得距离，得出最小的欧几里得距离

def euclidean_distance(p, q):  #求欧几里得距离
    p = np.array(p).flatten()
    q = np.array(q).flatten()
    return np.sqrt(np.sum(np.power((p-q) ,2)))

minDist = np.finfo("float").max
Q = np.dot(W.T, D - mean_X)
for i in xrange (len(projections)):
    dist = euclidean_distance( projections[i], Q)
    if dist < minDist :
    minDist = dist

如果要识别的图像是样本图像之一，那么求得的最小的欧几里得距离对应的样本图像与要识别的图像是同一个人。若要识别的图像非样本之一或根本不是人脸，我们就需要有一个阀值与求得的最小的欧几里得距离作比较，若在阀值之外，则可以判断要识别的图像不在样本中。关于如何设置阀值，本文不再讨论。

人脸识别的方法有很多种，基于特征脸的识别只是其中一种。但要实现一个可用的人脸识别，还有很多问题要考虑。另外PCA本身对某些特定情况的原始数据集也存在一些缺点。
至此，关于PCA，将不再深入探讨。
在PCA的学习过程中，深感一种技术应用的背后，必有惊艳的数学原理支撑，体会了一把数学之美:)，但因本人数学水平有限（后悔大学时没好好学），对以上理解必存在错漏和不详之处，所以也是很欢迎你的批评指正。

维基百科：特征脸
compute pca with this useful trick
Guide to face recognition with Python
FACE RECOGNITION HOMEPAGE
人脸图像数据库下载