PCA人臉識別學習及C語言實現
http://blog.csdn.net/jinshengtao/article/details/18599165
人臉識別主要方法：
.Eigenfaces，PCA（Principal Component Analysis），Turk and Pentland，1991
.Fisherfaces，LDA（Linear Discriminant Analysis），Belhumeur, Hespanha and Kriegman，1997
.LBPH，Local Binary Pattern Histograms，Ahonen, Hadid and Pietikäinen，2004
本文的目的，是結合人臉識別體驗一把PCA，體會其內涵：降維。另外文獻說，PCA的識別效果一般比神經網路ANN好。本文有20張人臉用於訓練，10張人臉用於測試。
訓練樣本和測試樣本來自：http://cswww.essex.ac.uk/mv/allfaces/faces94.zip1.PCA人臉識別方法
將PCA方法用於人臉識別，其實是假設所有的人臉都處於一個低維線性空間，而且不同的人臉在這個空間中具有可分性。其具體做法是由高維 影象空間經PCA變換後得到一組新的正交基，對這些正交基做一定的取捨，保留其中的一部分生成低維的人臉空間，也即是人臉的特徵子空間。PCA人臉識別演算法步驟包括：
a.人臉影象預處理【我沒做，人臉大小都是高200，寬180】
b.讀入人臉庫，訓練形成特徵子空間 【特徵值、特徵向量的求法，採用我上一篇文章的QR演算法】
c.把訓練影象和測試影象投影到上一步驟中的特徵子空間上 【矩陣相乘】
d.選擇一定的距離函式進行判別 【歐氏距離，挑最小的匹配】
2.PCA人臉識別流程
a.讀入人臉庫，讀入每一個二維的人臉影象並轉化為一維的向量，每個人選定一定數量的人臉照片構成訓練集【共20張】，則訓練集是一個36000*20的矩陣。測試集共10張影象，每次選一張，則測試集是一個36000*1的矩陣。
樣本集：
測試集：
程式碼：
- void load_data(double *T,IplImage *src,int k)
- {
- int i,j;
- //一副影象壓縮成一維的，存在T的一列裡
- for (i=0;i<IMG_HEIGHT;i++)
- {
- for (j=0;j<IMG_WIDTH;j++)
- {
- T[(i*IMG_WIDTH+j)*TRAIN_NUM+k-1]= (double)(unsigned char)src->imageData[i*IMG_WIDTH+j];
- }
- }
- }
b.計算 PCA變換的生成矩陣Q。首先計算訓練集的協方差矩陣X，其中x1,x2,...,xn為第i副影象的描述，即xi為一個36000*1的列向量。
由於這個矩陣太大36000*36000，求特徵值和特徵向量比較坑，所以改為求 P=XTX 的特徵向量和特徵值，且有如下性質：
設e是矩陣P的特徵值λ對應的特徵向量，則有：
這裡，X*e也是矩陣Q的特徵值λ對應的特徵向量，可以如此變換。
程式碼：
- void calc_mean(double *T,double *m)
- {
- int i,j;
- double temp;
- for (i=0;i<IMG_WIDTH*IMG_HEIGHT;i++)
- {
- temp=0;
- for (j=0;j<TRAIN_NUM;j++)
- {
- temp = temp + T[i*TRAIN_NUM+j];
- }
- m[i] = temp/TRAIN_NUM;
- }
- }
- void calc_covariance_matrix(double *T,double *L,double *m)
- {
- int i,j,k;
- double *T1;
- //T = T -m
- for (i=0;i<IMG_WIDTH*IMG_HEIGHT;i++)
- {
- for (j=0;j<TRAIN_NUM;j++)
- {
- T[i*TRAIN_NUM+j] = T[i*TRAIN_NUM+j] - m[i];
- }
- }
- T1 = (double *)malloc(sizeof(double)*IMG_HEIGHT*IMG_WIDTH*TRAIN_NUM);
- //L = T' * T
- matrix_reverse(T,T1,IMG_WIDTH*IMG_HEIGHT,TRAIN_NUM);
- matrix_mutil(L,T1,T,TRAIN_NUM,IMG_HEIGHT*IMG_WIDTH,TRAIN_NUM);
- free(T1);
- }
c.計算生成矩陣P的特徵值和特徵向量，並挑選合適的特徵值和特徵向量，構造特徵子空間變化矩陣。這裡P是實對稱矩陣，可以採用上一篇的方法，先進行Household變換將P變成三對角矩陣，然後使用QR迭代演算法求解特徵值和特徵向量，迭代次數60，誤差eps=0.000001，程式碼：
- void cstrq(double a[],int n,double q[],double b[],double c[])
- {
- int i,j,k,u,v;
- double h,f,g,h2;
- for (i=0; i<=n-1; i++)
- for (j=0; j<=n-1; j++)
- { u=i*n+j; q[u]=a[u];}
- for (i=n-1; i>=1; i--)
- { h=0.0;
- if (i>1)
- for (k=0; k<=i-1; k++)
- { u=i*n+k; h=h+q[u]*q[u];}
- if (h+1.0==1.0)
- { c[i]=0.0;
- if (i==1) c[i]=q[i*n+i-1];
- b[i]=0.0;
- }
- else
- { c[i]=sqrt(h);
- u=i*n+i-1;
- if (q[u]>0.0) c[i]=-c[i];
- h=h-q[u]*c[i];
- q[u]=q[u]-c[i];
- f=0.0;
- for (j=0; j<=i-1; j++)
- { q[j*n+i]=q[i*n+j]/h;
- g=0.0;
- for (k=0; k<=j; k++)
- g=g+q[j*n+k]*q[i*n+k];
- if (j+1<=i-1)
- for (k=j+1; k<=i-1; k++)
- g=g+q[k*n+j]*q[i*n+k];
- c[j]=g/h;
- f=f+g*q[j*n+i];
- }
- h2=f/(h+h);
- for (j=0; j<=i-1; j++)
- { f=q[i*n+j];
- g=c[j]-h2*f;
- c[j]=g;
- for (k=0; k<=j; k++)
- { u=j*n+k;
- q[u]=q[u]-f*c[k]-g*q[i*n+k];
- }
- }
- b[i]=h;
- }
- }
- for (i=0; i<=n-2; i++) c[i]=c[i+1];
- c[n-1]=0.0;
- b[0]=0.0;
- for (i=0; i<=n-1; i++)
- { if ((b[i]!=0.0)&&(i-1>=0))
- for (j=0; j<=i-1; j++)
- { g=0.0;
- for (k=0; k<=i-1; k++)
- g=g+q[i*n+k]*q[k*n+j];
- for (k=0; k<=i-1; k++)
- { u=k*n+j;
- q[u]=q[u]-g*q[k*n+i];
- }
- }
- u=i*n+i;
- b[i]=q[u]; q[u]=1.0;
- if (i-1>=0)
- for (j=0; j<=i-1; j++)
- { q[i*n+j]=0.0; q[j*n+i]=0.0;}
- }