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
人臉圖像數據庫下載