Programming Computer Vision with Python （學習筆記三）

概要

原書對於PCA的講解只有一小節，一筆帶過的感受，但我發現PCA是一個很重要的基礎知識點，在機器機視覺、人臉識別以及一些高級圖像處理技術時都被常常用到，因此本人自行對PCA進行了更深刻的學習。

PCA是什麼

PCA（Principal Component Analysis，主成分分析或主元分析）是一種算法，PCA的結果是用盡量少的特徵數據來表達最多的原始圖像的本質結構特徵。即便它會丟失一部分原始圖像的特徵表達，但它仍然是頗有用的處理技巧，也很經常使用，特別在計算機視覺和人臉識別方面。

假設咱們有一個二維數據集，它在平面上的分佈以下圖：

若是咱們想要用一個一維向量來表達此數據集，就會丟失一部分此數據集的信息，但咱們的目標是讓求得的這個一維向量能夠儘量多地保留這個數據集的特徵信息，那麼這個求解過程就是PCA。

經過PCA咱們能夠找到若干個1維向量，如圖：

直觀上看出，向量u1是數據集變化的主方向，而u2是次方向，u1比u2保留了更多的數據集的結構特徵，因此咱們選擇u1做爲主成分，並把原數據集投影到u1上就能夠得出對原數據集的一維近似重構：

以上只是一種直觀的示例，把二維數據降爲用1維來表示，固然，PCA一般是應用在高維數據集上。

PCA解決什麼問題

假設咱們有10張100 × 100像素的灰度人臉圖，咱們目標是要計算這10張圖的主成分來做爲人臉特徵，這樣就能夠基於這個‘特徵臉’進行人臉匹配和識別。但即便一個100 × 100像素的灰度圖像維度就達到10,000維，10張圖像的線性表示能夠達到100,000維，如此高維的數據帶來幾個問題：

• 對高維數據集進行分析處理的計算量是巨大的，消耗資源太大，時間太長

• 高維數據包含了大量冗餘和噪聲數據，會下降圖像識別率

因此一般對於高維數據集，首先須要對其進行降維運算，以低維向量表達原數據集最多最主要的結構特徵。從而將高維圖像識別問題轉化爲低維特徵向量的識別問題，大大下降了計算複雜度，同時也減小了冗餘信息所形成的識別偏差。PCA其實就是最經常使用的一種降維算法。

PAC也可用於高維數據壓縮、高維數據可視化（轉二維或三維後就能夠畫圖顯示）等方面，也是其它不少圖像處理算法的預處理步驟。

PCA的計算

關於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計算特徵臉（人臉識別中的一步），它應用在多張圖片上面。爲直觀起見，我用另一個例子——對單張黑白圖像進行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函數索引