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

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

1. 在计算协方差和特征向量的方法上，书上使用的是一种被做者称为compact trick的技巧，以及奇异值分解（SVD），这些都是什么东西呢？

2. 如何把PCA运用在多张图片上？

因此，咱们须要进一步的了解，同时，为示例对多张图片进行PCA，我选了一个跟书类似但更有趣的例子来作——人脸识别。

人脸识别与特征脸

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

特征脸计算步骤

1. 准备一组（假设10张）具备相同分辨率（假设:100 × 100）的人脸图像，把每张图像打平（numpy.flatten）成一个向量，即全部像素按行串联起来, 每一个图像被看成是10000维空间的一个点。再把全部打平的图像存储在 10000 × 10 矩阵X中，矩阵的每一列就是一张图片，每一个维为一行，共10000维

2. 对X的每一维（行）求均值，获得一个10000 × 1的向量，称为平均脸（由于把它变换为二维图像看起来像人脸），而后将X减去平均脸，即零均值化

3. 计算X的协方差矩阵C=XX^(X^表示X的转置)

4. 计算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(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、L一、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
人脸图像数据库下载